INTRODUCTION

El objetivo de este trabajo es realizar un análisis multivariante que nos permita indagar en la naturaleza de dos variables descritas por Brynjolfsson & McAfee (2014) en The Second Machine Age, como son el spread y el bounty. En concreto el spread creemos que será más sencillo de encontrar mediante un anáisis factorial que el bounty, ya que este último tiene una descripción algo más confusa. Podemos utilizar la descripción de esta sinopsis para definir ambas, pues resume muy bien ambos conceptos:

‘What the authors call “bounty” is their attempt to measure the benefits of new technology in ways reaching beyond such measures as GDP, which they say is inadequate. They use “spread” as a shorthand way to describe the increasing inequality that is also resulting from widespread new technology.’

Además, una vez hayamos conseguido definir bien la variable latente spread, podremos ver cómo están posicionados los países del mundo respecto a la misma, y cuáles están desmarcándose del resto positiva o negativamente debido al progreso tecnológico y a una creciente desigualdad en la repartición de los beneficios del mismo.

Para ello utilizaremos datos de distintos datasets, entre los cuales están UN, WHO, WB. En todos estos datasets existe una gran ausencia de información, una infinidad de valores ausentes y apenas son capaces de ponerse de acuerdo para utilizar los mismos nombres para los países. Todo esto dificultará la obtención del factor óptimo que nos de el spread (o la posición dentro de esa escala de spread de cada uno de los países).

Para tratar de hallar este spread, se utilizarán variables económicas (las que se pueda, ya que muchas que además son interesantes como el GINI ratio tienen muchos valores ausentes – una pena porque esta variable apunta en la misma dirección de los autores, hacia la desigualdad, nos da una idea de cómo es la desigualdad dentro de un mismo pais, lo cual nos podría ayudar para sacar tanto el spread como el bounty), de salud (acceso a sanidad y al agua, mortalidad infantil…), de igualdad (porcentaje de mujeres en el labour force…), etc. Buscamos tener variables que no midan sólo el potencial económico de un país, ya que como bien dicen los autores este potencial económico medido en términos de GDP es insuficiente a la hora de explicar el progresivo spread que se está generando entre unos países y otros en esta Second Machine Age. Habría sido muy positivo incluir variables que reflejen progreso tecnológico, venta de patentes etc; de hecho estas variables fueron descargadas del World Bank, pero tienen una infinidad de valores ausentes por lo que serán eliminadas al menos en su mayoría como puede verse en el código más adelante.

Otras medidas generales de progreso que investigaremos serán lo ecológico que es cada país en el uso de las energías, ya que como bien explican los autores, este “spread” debido al progreso tecnológico tendrá que ser sostenible para poder considerarse como tal, pues en el largo plazo no se mantendría de otra forma. El uso de internet y de celulares, la proporción de la economía dedicada a la agricultura, etc serán otras variables que se explorarán. Primero nos descargaremos todas las variables que nos encajen dentro de esta definición que dan Brynjolfsson & McAfee, bien sea para medir el spread o el bounty (aunque este ya decimos que es más difícil de medir, más aun sin poder contar con variables como el GINI ratio, que nos da la desigualdad de un país, que podría entenderse como la inversa de los beneficios del desarrollo tecnológico, tal y como lo describen ellos, ya que el bounty se encarga de explicar los beneficios para la población en general. De todas formas, creemos de partida que esta variable sólo podría entenderse como tal si pudiéramos hacer un análisis a lo largo del tiempo del factor que identifiquemos con la misma, ya que de otra forma no existe manera de constatar que efectivamente significa mejora de la calidad de vida debido al progreso tecnológico).

Después de descargar las variables y hacer un tratamiento exhaustivo de las mismas (por ejemplo armonizando los nombres y los códigos de los países en las diferentes bases de datos que nos bajemos), eliminaremos inteligentemente algunas observaciones; primero eliminamos las columnas que contienen demasiados valores ausentes, y después, sobre el data frame resultante, limpiamos los países que tienen información de pocas variables. De esta forma nos quedaremos con una matriz más completa, con un número inferior de valores ausentes, por lo que al utilizar la librería mice para la imputación de los mismos el random forest que utilizaremos para imputar variables tendrá suficientes observaciones cubiertas con respecto a las que tiene que rellenar. Al hacerlo de esta manera y no imputar directamente todo el dataset quedándonos con todas las variables que nos apetezca por muchos valores que tengan, reducimos la posibilidad de que se introduzca un sesgo en la muestra por haber utilizado pocas observaciones completas para rellenar muchos valores ausentes, en cuyo caso independientemente del modelo que uses te estás inventando los datos, con lo que el análisis completo quedaría invalidado.

Una vez tengamos el dataset completo con las variables y observaciones que queremos, procederemos a hacer un pequeño análisis descriptivo de las mismas; dibujaremos algunos gráficos, principalmente, para visualizar las distribuciones de variables, así como posibles outliers. También analizaremos la muestra de datos completa en su conjunto, viendo cómo se relacionan todas las variables con las demás. Por último, se realizará primero un PCA para ver cuánto porcentaje de la varianza podría coger con los primeros componentes, lo cual nos dará una idea también de qué tal seleccionadas están las variables. Posteriormente, trataremos de encontrar las variables spread y bounty realizando un Factor Analysis (toda la explicación respecto al número de factores que escoger etc se explicará más adelante). El modelo tiene la forma siguiente:

\(Variable_i = \mu_{i} + L_{i, j} \cdot Factor_j + \epsilon_{i}\)

\(i = 1, 2, ..., n\)

\(j = 1, 2, ..., m\)

n = número de observaciones

m = número de factor.

Finalmente se dibujarán algunos gráficos del mundo con los factores extraídos del análisis, y veremos si cuadra intuitivamente con la idea que tenemos de spread y bounty a partir de lo leído en The Second Machine Age

El trabajo está estructurado de la siguiente forma:

  1. Recogida, tratamiento y selección de variables para el análisis.

  2. Análisis descriptivo.

2.1. Análisis descriptivo individual.

2.2. Análisis descriptivo conjunto.

2.3. Inferencia.

  1. PCA & Factor Analysis: En la búsqueda del spread y el bounty.

3.1. PCA

3.2. FA

  1. Visualización de resultados,

1. RECOGIDA, TRATAMIENTO Y SELECCIÓN DE VARIABLES PARA EL ANÁLISIS.

En esta parte nos dedicaremos a leer y limpiar las variables. La mayoría de esta parte es código, ya que no hay mucho que comentar, pues no se centra en el análisis. En algunos casos cogeremos las variables para el último año, mientras que en otros casos tendremos que usar la media de los últimos años, o datos de 2016. Todo esto se verá explícitamente en el código. El procedimiento general es el siguiente: leemos la variable, utilizamos dplyr para dejarla en el mismo formato que a las demás, cambiamos nombres de columnas, seleccionamos solo las columnas que nos queremos quedar, identificamos de qué forma vienen los NAs para codificar valores como “…” como NAs, y hacemos merge para juntar cada variable a todas las demás.

Como hemos seleccionado sólo unas pocas de las variables totales consideradas, sólo las describiremos e indicaremos de dónde nos las bajamos de estas variables, para no alargar excesivamente sin sentido el trabajo.

Este dataset de aquí abajo está bajado de aquí, y contiene datos recopilados por las naciones unidas de todos los países del mundo. Nos quedaremos de ese dataset sólo con las variables que nos resulten de especial interés. Sin embargo, debemos decir que es un dataset mucho más completo que los encontrados en World Health Organization y en el World Bank; para las mismas variables tienen en la mayoría de los casos muchas más observaciones, no contienen tantos datos ausentes y es además un mapa más completo del mundo, al contener muchos más países.

Como los nombres de países no son homogéneos, trataremos de armonizar los nombres de estos en los distintos datasets de la siguiente forma:

Country_Name_correct<- c("Antigua and Barbuda"="Antigua", 
                         "Bahamas, The"="Bahamas", 
                         "Brunei Darussalam"="Brunei",
                         "Cabo Verde"="Cape Verde", 
                         "Congo, Dem. Rep."="Democratic Republic of the Congo", 
                         "Congo, Rep."="Republic of Congo", 
                         "Cote d'Ivoire"="Ivory Coast", 
                         "Egypt, Arab Rep."="Egypt", 
                         "Faeroe Islands"="Faroe Islands", 
                         "Gambia, The"="Gambia", 
                         "Iran, Islamic Rep."="Iran", 
                         "Korea, Dem. Rep."="North Korea", 
                         "Korea, Rep."="South Korea", 
                         "Korea, Republic"="South Korea", 
                         "Republic of Korea"="South Korea",
                         "Kyrgyz Republic"="Kyrgyzstan", 
                         "Lao PDR"="Laos", 
                         "Macedonia, FYR"="Macedonia, Federal States of", 
                         "Micronesia, Fed. Sts."="Micronesia, Federal States of", 
                         "Russian Federation"="Russia", 
                         "Slovak Republic"="Slovakia", 
                         "St. Lucia"="Saint Lucia", 
                         "St. Martin (French part)"="Saint Martin", 
                         "St. Vincent and the Grenadines"="Saint Vincent",
                         "Saint Vincent and the Grenadines"="Saint Vincent", 
                         "Syrian Arab Republic"="Syria", 
                         "Trinidad and Tobago"="Trinidad", 
                         "United Kingdom"="UK",
                         "United States"="USA", 
                         "United States of America" = "USA", 
                         "Venezuela, RB"="Venezuela", 
                         "Virgin Islands (U.S.)"="Virgin Islands",
                         "United States Virgin Islands"="US Virgin Islands", 
                         "Yemen, Rep."="Yemen",
                         "Sint Maarten (Dutch part)"="Sint Maarten",
                         "West Bank and Gaza"="Palestine",
                         "The former Yugoslav Republic of Macedonia" = "Macedonia",
                         "Democratic People's Republic of Korea"="North Korea",
                         "Bolivia (Plurinational State of)"="Bolivia",
                         "Lao People's Democratic Republic"="Laos", 
                         "Micronesia (Federated States of)"="Micronesia, Federal States of", 
                         "Venezuela (Bolivarian Republic of)"="Venezuela", 
                         "State of Palestine" = "Palestine", 
                         "China, Hong Kong SAR"="Hong Kong", 
                         "China, Macao SAR" = "Macao SAR, China", 
                         "Iran (Islamic Republic of)"="Iran", 
                         "Korea, Dem. People's Rep."="North Korea", 
                         "United Republic of Tanzania"="Tanzania", 
                         "Republic of Moldova"="Moldova", 
                         "Viet Nam" = "Vietnam", 
                         "Virgin Islands"="US Virgin Islands", 
                         "Hong Kong SAR"="Hong Kong", 
                         "Micronesia" = "Micronesia, Federal States of", 
                         "Czech Republic" = "Czechia", 
                         "Macao" = "Macao SAR, China",
                         "Congo" = "Republic of Congo", 
                         "Saint Kitts and Nevis"="St. Kitts and Nevis",
                         "Swaziland"="Eswatini")
                         
                         
correct_names <- function(df, correct_list = Country_Name_correct) {
  
  df2 <- df
  
  for (c in names(Country_Name_correct) ) {
    
    df2[df2$country==c,"country"] <-  Country_Name_correct[c]
  
  }
  return(df2)
}



add_code_absent <- function(df, countries = c("Eswatini", "Micronesia"), codes = c("SWZ", "FSM")) { #estos dos países nos dan fallos al meterles el código, así que usaremos la información de: https://www.iso.org/obp/ui/#iso:code:3166:SZ
  
  for(i in 1:length(countries)) {
    
    df[df$country == countries[i], "country_code"] <- codes[i]
    
  }
  
  return(df)
}


if(!require(countrycode)) {
  
  install.packages("countrycode")
  
}
## Loading required package: countrycode
library(countrycode)
country_profile_vars <- read.csv("country_profile_variables.csv", sep = ",", fileEncoding = "utf-8", stringsAsFactors = F)

country_profile_vars <- correct_names(country_profile_vars)

country_profile_vars$country_code <- countrycode(country_profile_vars$country, "country.name", "iso3c")
## Warning in countrycode(country_profile_vars$country, "country.name", "iso3c"): Some values were not matched unambiguously: Channel Islands, Eswatini
country_profile_vars <- add_code_absent(country_profile_vars)

country_profile_vars <- country_profile_vars[!is.na(country_profile_vars$country_code), ]

countries_and_codes_total <- country_profile_vars[!duplicated(country_profile_vars$country_code) , c("country", "country_code")]
internet_usage <- country_profile_vars[ , c("country", "country_code", "Individuals.using.the.Internet..per.100.inhabitants.")] 


names(internet_usage) <- c("country", "country_code","internet_usage")

internet_usage[internet_usage$internet_usage == -99, "internet_usage"] <- NA

internet_usage <- correct_names(internet_usage)

Echando un vistazo a la variable, salta a la vista que algunos valores están mal introducidos, ya que es imposible que, haciéndolo relastivo a número de personas, en Brasil se acceda más a Internet que en España, pues es bien conocido que es un país con una buena parte de la población viviendo en la pobreza. Por este motivo muy probablemente acabe eliminada.

mobile_cellular_subscriptions <- country_profile_vars[ , c("country" ,"country_code", "Mobile.cellular.subscriptions..per.100.inhabitants.")]

names(mobile_cellular_subscriptions) <-  c("country", "country_code","cell_subscr")

mobile_cellular_subscriptions$cell_subscr[mobile_cellular_subscriptions$cell_subscr == -99 | mobile_cellular_subscriptions$cell_subscr == "..."] <- NA

#mobile_cellular_subscriptions[mobile_cellular_subscriptions$country %in% inf_mort_rate$country, ]    

mobile_cellular_subscriptions <- correct_names(mobile_cellular_subscriptions)

internet_usage$country <- NULL

mobile_cellular_subscriptions$country <- NULL

mydata <- merge(internet_usage, mobile_cellular_subscriptions, on = "country_code", all = T)
gdp <- data.frame(country_profile_vars$country, country_profile_vars$country_code ,country_profile_vars$GDP..Gross.domestic.product..million.current.US..)

names(gdp) <- c("country","country_code", "gdp")

gdp$gdp[gdp$gdp == -99 | gdp$gdp == "..."] <- NA

gdp <- gdp[!is.na(gdp$country), ]

gdp$country <- as.character(gdp$country)

gdp <- correct_names(gdp)

gdp$country <- NULL

mydata <- merge(mydata, gdp, on = "country_code", all = T)
library(readr)

worldbankdata <- read_csv("worldbank1.csv")
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `GDP per capita (Current US$)` = col_character(),
##   `Labor force, female (% of total labor force)` = col_character(),
##   `Inflation, consumer prices (annual %)` = col_character(),
##   `GNI per capita, Atlas method (current US$)` = col_character(),
##   `Maternity leave (days paid)` = col_character(),
##   `Ratio of female to male labor force participation rate (%) (modeled ILO estimate)` = col_character(),
##   `Contributing family workers, female (% of female employment) (modeled ILO estimate)` = col_character()
## )
worldbankdata <- worldbankdata[ , c("Country Name", "GDP per capita (Current US$)", "Labor force, female (% of total labor force)",  "Inflation, consumer prices (annual %)", "GNI per capita, Atlas method (current US$)", "Maternity leave (days paid)", "Ratio of female to male labor force participation rate (%) (modeled ILO estimate)", "Contributing family workers, female (% of female employment) (modeled ILO estimate)")]

names(worldbankdata) <- c("country", "gdp_per_cap", "lab_force_fem", "inflation", "gni_per_cap", "matern_leave", "fem_to_male_part_rate", "contr_fam_workers")

worldbankdata[worldbankdata == ".."] <-  NA

worldbankdata <- correct_names(worldbankdata)

worldbankdata$country_code <- countrycode(worldbankdata$country, "country.name", "iso3c")
## Warning in countrycode(worldbankdata$country, "country.name", "iso3c"): Some values were not matched unambiguously: Caribbean small states, Channel Islands, Eswatini, Euro area, Kosovo, North America, Saint Martin, Small states, South Asia, World
worldbankdata <- add_code_absent(worldbankdata)

worldbankdata <- worldbankdata[!is.na(worldbankdata$country_code), ]

worldbankdata <- data.frame(worldbankdata)

worldbankdata[ , 2:8] <- apply(worldbankdata[ , 2:8], 2, as.numeric)

worldbankdata$country <- NULL

mydata <- merge(mydata, worldbankdata, on = "country_code", all = T)
balance_of_payments <- country_profile_vars[ , c("country", "country_code","Balance.of.payments..current.account..million.US..")]

names(balance_of_payments) <- c("country","country_code", "BoP")

balance_of_payments$BoP[balance_of_payments$BoP == -99 | balance_of_payments$BoP == "..."] <- NA

balance_of_payments <- correct_names(balance_of_payments)

balance_of_payments$country <- NULL

mydata <- merge(mydata, balance_of_payments, on = "country_code", all = T)
urban_pop <- country_profile_vars[, c("country", "country_code","Urban.population....of.total.population.")]

names(urban_pop) <- c("country", "country_code","urban_p")

urban_pop$urban_p[urban_pop$urban_p == -99 | urban_pop$urban_p == ".."] <- NA

urban_pop <- correct_names(urban_pop)

urban_pop$country <- NULL

mydata <- merge(mydata, urban_pop, on = "country_code", all = T)
refugees <- country_profile_vars[ , c("country", "country_code","Refugees.and.others.of.concern.to.UNHCR..in.thousands." )]

names(refugees) <- c("country", "country_code","refugees")

refugees$refugees[refugees$refugees == -99 | refugees$refugees  == ".."] <- NA

refugees <- correct_names(refugees)

refugees[refugees == "~0.0"] <- NA

refugees$country <- NULL

mydata <- merge(mydata,refugees, on = "country_code", all = T)
gini_index <- read_csv("gini_index.csv")

names(gini_index)
##  [1] "Series Name"   "Series Code"   "Country Name"  "Country Code" 
##  [5] "1990 [YR1990]" "2000 [YR2000]" "2008 [YR2008]" "2009 [YR2009]"
##  [9] "2010 [YR2010]" "2011 [YR2011]" "2012 [YR2012]" "2013 [YR2013]"
## [13] "2014 [YR2014]" "2015 [YR2015]" "2016 [YR2016]" "2017 [YR2017]"
country_and_code <- gini_index[ , c("Country Name", "Country Code")]

names(country_and_code) <- c("country", "country_code")

gini_index<- gini_index[ ,9:ncol(gini_index)]

gini_index <- apply(gini_index, 2, as.numeric)

gini_index$mean <- rowMeans(gini_index, na.rm = T)

gini_index <- data.frame(cbind(country_and_code, gini_index$mean))

names(gini_index) <- c("country", "country_code", "gini_index")

gini_index <- correct_names(gini_index)

gini_index$country <- NULL

mydata <- merge(mydata, gini_index, on = "country_code", all = T)
health <- country_profile_vars[, c("country", "country_code","Health..Total.expenditure....of.GDP.")]

names(health) <- c("country", "country_code","health")

health$health[health$health == -99 | health$health  == ".."] <- NA

health <- correct_names(health)

health$country <- NULL

mydata <- merge(mydata,health, on = "country_code", all = T)
education <- country_profile_vars[ , c("country", "country_code","Education..Government.expenditure....of.GDP.")]

names(education) <- c("country", "country_code","education")

education$education[education$education == -99 | education$education  == ".." | education$education  == "..."] <- NA

education <- correct_names(education)

education$country <- NULL

mydata <- merge(mydata,education, on = "country_code", all = T)
co2 <- country_profile_vars[ , c("country", "country_code","CO2.emission.estimates..million.tons.tons.per.capita.")]

names(co2) <- c("country", "country_code","co2")

co2$co2[co2$co2 == -99 | co2$co2  == ".."] <- NA

co2 <- correct_names(co2)

co2$country <- NULL

mydata <- merge(mydata,co2, on = "country_code", all = T)
productivity <- read_csv("productivity_wb.csv") #as measured by GDP per person employed. 
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `GDP per person employed (constant 2011 PPP $)` = col_character()
## )
productivity <- productivity[ , c(3, 4, 5)]

names(productivity) <- c("country", "country_code", "productivity")

productivity$productivity[productivity$productivity == -99 | productivity$productivity == ".."] <- NA

productivity <- correct_names(productivity)

productivity <- add_code_absent(productivity)

productivity$country <- NULL

mydata <- merge(mydata, productivity, on = "country_code", all = T)
innovation_and_fin_health <- read_csv("innovation_and_debt.csv")
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Scientific and technical journal articles` = col_character(),
##   `Patent applications, nonresidents` = col_character(),
##   `Patent applications, residents [IP.PAT.RESD]` = col_character(),
##   `Trademark applications, total` = col_character(),
##   `High-technology exports (current US$)` = col_character(),
##   `Charges for the use of intellectual property, receipts (BoP, current US$)` = col_character(),
##   `Depth of credit information index (0=low to 8=high)` = col_character(),
##   `Depth of the food deficit (kilocalories per person per day)` = col_character(),
##   `IMF charges (INT, current US$)` = col_character(),
##   `Interest forgiven (current US$)` = col_character()
## )
head(innovation_and_fin_health)
## # A tibble: 6 x 14
##    Time `Time Code` `Country Name` `Country Code` `Scientific and~
##   <int> <chr>       <chr>          <chr>          <chr>           
## 1  2016 YR2016      Afghanistan    AFG            80.4            
## 2  2016 YR2016      Albania        ALB            191.4           
## 3  2016 YR2016      Algeria        DZA            4447.1          
## 4  2016 YR2016      American Samoa ASM            ..              
## 5  2016 YR2016      Andorra        AND            7.9             
## 6  2016 YR2016      Angola         AGO            39.1            
## # ... with 9 more variables: `Patent applications, nonresidents` <chr>,
## #   `Patent applications, residents [IP.PAT.RESD]` <chr>, `Trademark
## #   applications, total` <chr>, `High-technology exports (current
## #   US$)` <chr>, `Charges for the use of intellectual property, receipts
## #   (BoP, current US$)` <chr>, `Depth of credit information index (0=low
## #   to 8=high)` <chr>, `Depth of the food deficit (kilocalories per person
## #   per day)` <chr>, `IMF charges (INT, current US$)` <chr>, `Interest
## #   forgiven (current US$)` <chr>
innovation_and_fin_health <- innovation_and_fin_health[ , 3:ncol(innovation_and_fin_health)]

names(innovation_and_fin_health) <- c("country", "country_code", "sci_jour", "patent_nr", "patent_r", "trademark app", "hightech_exp", "int_prop", "credit_inf", "food_deficit", "imf_charges", "int_forgiven")

innovation_and_fin_health[innovation_and_fin_health == ".." | innovation_and_fin_health == -99] <- NA

innovation_and_fin_health <- correct_names(innovation_and_fin_health)

innovation_and_fin_health <- add_code_absent(innovation_and_fin_health)

innovation_and_fin_health$country <- NULL

mydata <- merge(mydata, innovation_and_fin_health, on = "country_code", all = T)
renewable_energy <- read_csv("renewable_energy.csv")
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Renewable energy share of TFEC (%)` = col_character(),
##   `Renewable electricity share of total electricity output (%)` = col_character()
## )
renewable_energy <- renewable_energy[ , 3:ncol(renewable_energy)]

names(renewable_energy) <- c("country", "country_code", "renew_en", "renew_electr")

renewable_energy <- correct_names(renewable_energy)

renewable_energy <- add_code_absent(renewable_energy)

renewable_energy$country <- NULL

mydata <- merge(mydata, renewable_energy, on = "country_code", all = T)
library(WHO)
library(rgho)

library(plyr)
library(dplyr)
require(tidyr)

infant_mortality_rate <- data.frame(get_data("MDG_0000000001"))

inf_mort_rate <-  infant_mortality_rate %>% 
  
  group_by(country) %>% 
  
    group_by(year, add = T) %>% 
  
      filter(year >= 2010) %>% 
  
        select(country, value) %>%
  
          filter(!is.na(country)) %>%
  
            spread(year, value) 
  
inf_mort_rate$inf_mort_rate <- inf_mort_rate$`2017`

inf_mort_rate <- inf_mort_rate[ , c("country", "inf_mort_rate")]


inf_mort_rate <- correct_names(inf_mort_rate)

inf_mort_rate$country_code <- countrycode(inf_mort_rate$country, "country.name", "iso3c")

inf_mort_rate <- add_code_absent(inf_mort_rate)
sanitation_services <- read.csv("sanitation_services.csv")

names(sanitation_services) <- c("country", "sanitation_services")

sanitation_services <- sanitation_services[!is.na(sanitation_services$country), ]

sanitation_services$country <- as.character(sanitation_services$country)

sanitation_services <- correct_names(sanitation_services)

sanitation_services$country_code <- countrycode(sanitation_services$country, "country.name", "iso3c")

sanitation_services <- add_code_absent(sanitation_services)

inf_mort_rate$country <- NULL

sanitation_services$country <- NULL

who_df <- merge(inf_mort_rate, sanitation_services, on = "country_code", all = T)
deaths_for_household_pollution <- read.csv("deaths_for_household_pollution.csv", sep = ",", fileEncoding = "utf-8", stringsAsFactors = F)

deaths_for_household_pollution$Value <- as.numeric(gsub(' .*', "", deaths_for_household_pollution$Value))

deaths_for_household_pollution <- deaths_for_household_pollution[deaths_for_household_pollution$Cause == "Total", ]

names(deaths_for_household_pollution) <- c("country", "cause", "deaths_household_pollution")

deaths_for_household_pollution$cause <- NULL

deaths_for_household_pollution <- deaths_for_household_pollution[!is.na(deaths_for_household_pollution$country), ]

deaths_for_household_pollution <- correct_names(deaths_for_household_pollution)

deaths_for_household_pollution$country_code <- countrycode(deaths_for_household_pollution$country, "country.name", "iso3c")

deaths_for_household_pollution <- add_code_absent(deaths_for_household_pollution)

deaths_for_household_pollution$country <- NULL

who_df <- merge(who_df, deaths_for_household_pollution, on = "country_code", all = T)
water <- read.csv("water.csv", stringsAsFactors = F)

names(water) <- c("country", "water")

water <- water[!is.na(water$country), ]

water <- correct_names(water)

water$country_code <- countrycode(water$country, "country.name", "iso3c")

water <- add_code_absent(water)

water$country <- NULL

who_df <- merge(who_df, water, on = "country_code", all = T)
require(tidyr)
library(plyr)
require(dplyr)

energy_wb <- read_csv("energy_wb.csv")
## Parsed with column specification:
## cols(
##   Time = col_character(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Access to electricity (% of total population) [1.1_ACCESS.ELECTRICITY.TOT]` = col_character(),
##   `Access to Clean Fuels and Technologies for cooking (% of total population) [2.1_ACCESS.CFT.TOT]` = col_character()
## )
energy_wb <- energy_wb[1:460, c(1, 3:6)]

names(energy_wb) <- c("year", "country", "country_code", "electricity", "clean_fuels")


library(data.table) ## v >= 1.9.6
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
energy_wb <- dcast(setDT(energy_wb), country ~ year, value.var = c("electricity", "clean_fuels"))

energy_wb[energy_wb == ".." | energy_wb == -99] <- NA


energy_new <- data.frame(energy_wb)

for(column in c("electricity", "clean_fuels")) {
  
  nombre_real <- paste0(column, "_2016")
  
  for (i in 1:nrow(energy_wb)) {
    
    if(is.na(energy_new[i, nombre_real])) {
      
      if(!is.na(energy_new[i, paste0(column, "_2015")])) {
        
        energy_new[i, nombre_real] <- energy_new[i, paste0(column, "_2015")]
        
      } else {
        
        next
        
      }
      
    }
    
  }
  
}

energy_wb <- correct_names(energy_wb)

energy_new$country_code <- countrycode(energy_new$country, "country.name", "iso3c")
## Warning in countrycode(energy_new$country, "country.name", "iso3c"): Some values were not matched unambiguously: Channel Islands, Kosovo, Netherlands Antilles
energy_new <- add_code_absent(energy_new)

energy_new <- energy_new[!is.na(energy_new$country_code), ]

energy_new$electricity_2015 <- NULL

energy_new$clean_fuels_2015 <- NULL

energy_new$country <- NULL

mydata <- merge(mydata, energy_new, on = "country_code", all = T)
mydata <- merge(mydata, who_df, on = "country_code", all = T)

#codes_and_countries <- mydata[ , c("country", "country_code")]
#mydata$country_code <- NULL

En este chunk de abajo se irán eliminando algunas variables que, después de haber realizado el análisis de factores completo y de haber visualizado suficientemente la relación multidimensional con el resto de variables, resultan no tener mucho que aportar a la construcción del modelo, y hacen que este no pueda captar la mayoría de la varianza total de la muestra de variables.

additional_un_vars <- country_profile_vars[ , c("country_code", "Fertility.rate..total..live.births.per.woman.", "Life.expectancy.at.birth..females.males..years.", "Urban.population.growth.rate..average.annual...", "International.trade..Balance..million.US..", "Employment..Agriculture....of.employed.", "Economy..Agriculture....of.GVA.", "Economy..Industry....of.GVA.")]

additional_un_vars[additional_un_vars== ".." | additional_un_vars == "..." | additional_un_vars == -99 | additional_un_vars == ".../..."] <- NA

names(additional_un_vars) <- c("country_code", "fert_rate", "life_exp", "urban_pop_growth", "int_tr_bal", "emp_agr", "eco_agr", "eco_ind")

life_exp <- c()

for (i in 1:nrow(additional_un_vars)) {
  
  if(!is.na(additional_un_vars[i, "life_exp"])) {
    
      fem <- as.numeric(gsub("/.*", "", additional_un_vars[i, "life_exp"]))
      
      mas <- as.numeric(gsub(".*/", "", additional_un_vars[i, "life_exp"]))
      
      media <- (fem + mas) / 2
      
      life_exp <- c(life_exp, media)
      
  } else {
    
    life_exp <- c(life_exp, NA)
    
  }
  
}

additional_un_vars$life_exp <- life_exp


mydata <- merge(mydata, additional_un_vars, on = "country_code", all = T)
vulner_empl <- read_csv("vulnerable_employment.csv")
## Parsed with column specification:
## cols(
##   Time = col_character(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Vulnerable employment, total (% of total employment) (modeled ILO estimate) [SL.EMP.VULN.ZS]` = col_character()
## )
vulner_empl <- vulner_empl[1:217 , 3:ncol(vulner_empl)]

names(vulner_empl) <- c("country", "country_code", "vulner_empl")

vulner_empl[vulner_empl == ".." | vulner_empl == -99] <- NA

vulner_empl$country <- NULL

mydata <- merge(mydata, vulner_empl, on = "country_code", all = T)

La variable domestic credit viene de World Bank y es una variable que nos ayuda a entender el nivel de confianza entre los agentes económicos de un país. Esperamos que cuanto mayor sea el nivel de spread, o cuanto mejor posicionados estén los países en la escala de spread, mayor será también la confianza, identificada en este caso como el credito doméstico al sector privado como porcentaje del GDP en el año 2017.

dom_credit <- read_csv("dom_credit.csv")
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Domestic credit to private sector (% of GDP) [FS.AST.PRVT.GD.ZS]` = col_character()
## )
dom_credit <- dom_credit[ , 3:5]

names(dom_credit) <- c("country", "country_code", "dom_credit")

dom_credit[dom_credit == "..." | dom_credit == ".." | dom_credit == -99] <- NA

dom_credit$country <- NULL

mydata <- merge(mydata, dom_credit, on = "country_code", all = T)

La variable account ownership es una variable que nos indica el porcentaje de personas de un país que tienen acceso, bien sea por tener una cuenta en una entidad financiera o bien por gozar de este servicio vía móvil, al sistema económico y financiero; cuántas personas forman parte de ese sistema dentro de un país. Intuitivamente, en los países que mejor estén aprovechando este progreso tecnológico habrá actividad económica, y esta variable nos dice cuánta gente participa en esa actividad económica; este es un asunto comentado por los autores de especial importancia a la hora de ir “beyond gdp”, pues para que el bounty sea sostenible, y puedas posicionarte como país arriba en la escala de spread, este bounty tiene que estar bien distribuido entre tu poblacón; de otra forma no existirá equilibrio social ni económico.

acc_own <- read_csv("account_ownership.csv")
## Parsed with column specification:
## cols(
##   Time = col_integer(),
##   `Time Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `Account ownership at a financial institution or with a mobile-money-service provider (% of population ages 15+) [FX.OWN.TOTL.ZS]` = col_character()
## )
acc_own <- acc_own[ , 3:5]

names(acc_own) <- c("country", "country_code", "acc_own")

acc_own[acc_own == "..." | acc_own == ".." | acc_own == -99] <- NA

acc_own$country <- NULL

mydata <- merge(mydata, acc_own, on = "country_code", all = T)
educ_staff_comp <- read_csv("educ_staff_comp.csv")
## Parsed with column specification:
## cols(
##   `Series Name` = col_character(),
##   `Series Code` = col_character(),
##   `Country Name` = col_character(),
##   `Country Code` = col_character(),
##   `2010 [YR2010]` = col_character(),
##   `2011 [YR2011]` = col_character(),
##   `2012 [YR2012]` = col_character(),
##   `2013 [YR2013]` = col_character(),
##   `2014 [YR2014]` = col_character(),
##   `2015 [YR2015]` = col_character(),
##   `2016 [YR2016]` = col_character(),
##   `2017 [YR2017]` = col_character()
## )
educ_staff_comp <- educ_staff_comp[ , 3:ncol(educ_staff_comp)]

educ_staff_comp$ed_comp <- rowMeans(apply(educ_staff_comp[ ,3:ncol(educ_staff_comp)], 2, as.numeric ), na.rm=T)
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción

## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
educ_staff_comp <- educ_staff_comp[ , c(1,2,11)]

names(educ_staff_comp) <- c("country", "country_code", "ed_comp")

educ_staff_comp[educ_staff_comp == "..." | educ_staff_comp == ".." | educ_staff_comp == -99] <- NA

educ_staff_comp <-  educ_staff_comp[!is.na(educ_staff_comp$country), ]

educ_staff_comp$country <- NULL

mydata <- merge(mydata, educ_staff_comp, on = "country_code", all = T)
prop_miss <- function(x){sum(is.na(x))/length(x)*100}
apply(mydata,2,prop_miss)
##               country_code             internet_usage 
##                   0.000000                   3.361345 
##                cell_subscr                        gdp 
##                  11.344538                  11.764706 
##                gdp_per_cap              lab_force_fem 
##                  21.008403                  21.848739 
##                  inflation                gni_per_cap 
##                  31.092437                  22.689076 
##               matern_leave      fem_to_male_part_rate 
##                  25.210084                  21.428571 
##          contr_fam_workers                        BoP 
##                  21.428571                  23.949580 
##                    urban_p                   refugees 
##                   3.361345                  31.092437 
##                 gini_index                     health 
##                  42.857143                  19.327731 
##                  education                        co2 
##                  36.974790                  11.344538 
##               productivity                   sci_jour 
##                  22.268908                  16.806723 
##                  patent_nr                   patent_r 
##                  51.680672                  54.621849 
##              trademark app               hightech_exp 
##                  42.857143                  43.697479 
##                   int_prop                 credit_inf 
##                  49.579832                  19.747899 
##               food_deficit                imf_charges 
##                  50.840336                  47.478992 
##               int_forgiven                   renew_en 
##                  47.478992                   2.521008 
##               renew_electr           electricity_2016 
##                   2.521008                  10.084034 
##           clean_fuels_2016              inf_mort_rate 
##                  20.588235                  17.647059 
##        sanitation_services deaths_household_pollution 
##                  18.487395                  23.529412 
##                      water                  fert_rate 
##                  18.487395                  10.504202 
##                   life_exp           urban_pop_growth 
##                  10.924370                   3.361345 
##                 int_tr_bal                    emp_agr 
##                  10.924370                  18.067227 
##                    eco_agr                    eco_ind 
##                  12.605042                  11.764706 
##                vulner_empl                 dom_credit 
##                  20.588235                  41.596639 
##                    acc_own                    ed_comp 
##                  39.495798                  44.537815
apply(mydata,1,prop_miss)
##   [1] 54.166667 14.583333 16.666667 70.833333  2.083333 47.916667 93.750000
##   [8] 18.750000  2.083333  2.083333 83.333333 37.500000  8.333333  6.250000
##  [15]  4.166667 14.583333  6.250000  6.250000 87.500000  8.333333  2.083333
##  [22]  2.083333 12.500000 18.750000  6.250000  2.083333 12.500000 56.250000
##  [29]  0.000000  0.000000 16.666667 18.750000 14.583333  6.250000 31.250000
##  [36]  8.333333  8.333333 89.583333  6.250000  6.250000 41.666667  8.333333
##  [43] 16.666667 20.833333 58.333333  2.083333 18.750000 12.500000  2.083333
##  [50] 31.250000 85.416667 70.833333  6.250000  6.250000  6.250000 16.666667
##  [57] 39.583333  6.250000  6.250000  4.166667  4.166667  6.250000 33.333333
##  [64] 81.250000  6.250000  8.333333 10.416667  6.250000 10.416667 81.250000
##  [71]  6.250000 75.000000 43.750000 14.583333  6.250000  2.083333  2.083333
##  [78] 79.166667 12.500000 79.166667 12.500000 10.416667 29.166667 10.416667
##  [85] 29.166667 62.500000  2.083333 81.250000 68.750000 10.416667 14.583333
##  [92] 31.250000  2.083333  6.250000 14.583333  6.250000  2.083333 83.333333
##  [99]  4.166667  6.250000  8.333333 18.750000  8.333333  8.333333  6.250000
## [106]  8.333333  4.166667 10.416667  0.000000  2.083333  0.000000  6.250000
## [113] 39.583333 45.833333  4.166667  4.166667 18.750000 10.416667 12.500000
## [120] 12.500000 31.250000 18.750000 62.500000  4.166667 12.500000  6.250000
## [127]  8.333333  6.250000 43.750000 91.666667  4.166667 56.250000  2.083333
## [134]  0.000000 14.583333  0.000000 47.916667 10.416667 10.416667  6.250000
## [141]  8.333333 10.416667  2.083333 79.166667  4.166667 12.500000 72.916667
## [148] 79.166667  4.166667  6.250000  2.083333 81.250000 12.500000 52.083333
## [155]  6.250000 14.583333 12.500000 77.083333  8.333333  8.333333  4.166667
## [162] 60.416667 14.583333 16.666667  2.083333  4.166667  2.083333  6.250000
## [169] 45.833333 20.833333  8.333333 47.916667 43.750000  8.333333  8.333333
## [176] 35.416667 52.083333 16.666667 93.750000  2.083333  2.083333  2.083333
## [183] 14.583333 18.750000  8.333333 12.500000 81.250000 14.583333  4.166667
## [190]  0.000000 43.750000 33.333333 83.333333  4.166667 27.083333 12.500000
## [197] 20.833333  6.250000 10.416667  6.250000 10.416667 70.833333 27.083333
## [204] 29.166667 70.833333 14.583333  8.333333  0.000000 12.500000 87.500000
## [211] 22.916667 16.666667 25.000000 14.583333  2.083333  6.250000 52.083333
## [218] 93.750000 10.416667  2.083333  6.250000 10.416667 10.416667 18.750000
## [225] 91.666667 18.750000 25.000000 70.833333 68.750000  4.166667 18.750000
## [232] 85.416667 14.583333 79.166667 18.750000  2.083333 12.500000  4.166667
get_noisy_countries <- function(df) {
  
  require(lava)
  
  if ("country" %ni% names(df)) {
    
    noisy_countries <- c()
  
    for(i in 1:nrow(df)) {
    
      if(sum(is.na(df[i, ])) / ncol(df) > 0.3) {
      
        noisy_countries <- c(noisy_countries, df[i, "country_code"])
      
    } else {
      
      next
      
    }
    
    
   }
    
  } else {
    
      noisy_countries <- c()
    
      for(i in 1:nrow(df)) {
    
        if(sum(is.na(df[i, ])) / ncol(df) > 0.3) {
      
          noisy_countries <- c(noisy_countries, df[i, "country"])
      
    } else {
      
      next
      
    }
    
   }
    
  }

    return(noisy_countries)
  
}    

noisy_c <- get_noisy_countries(mydata)
## Loading required package: lava
## lava version 1.6.3
## 
## Attaching package: 'lava'
## The following object is masked from 'package:dplyr':
## 
##     vars
noisy_c
##  [1] "ABW" "AIA" "AND" "ANT" "ASM" "ATG" "BES" "BMU" "CAF" "CHI" "CIV"
## [12] "COK" "CUB" "CUW" "CYM" "DMA" "ERI" "ESH" "FLK" "FRO" "FSM" "GIB"
## [23] "GLP" "GRL" "GUF" "GUM" "HKG" "IMN" "KIR" "KNA" "LBY" "LIE" "MAC"
## [34] "MAF" "MCO" "MHL" "MNP" "MSR" "MTQ" "MYT" "NCL" "NIU" "NRU" "PLW"
## [45] "PRI" "PRK" "PSE" "PYF" "REU" "SHN" "SMR" "SOM" "SPM" "SXM" "TCA"
## [56] "TKL" "TUV" "TWN" "VAT" "VGB" "VIR" "WLF" "XKX"

Vamos a probar eliminando primero las variables que tengan muchos valores ausentes, y después volvemos a ver cuántos países ruidosos tenemos.

remove_noisy_vars <- function(df) {
  
  cols_pos <- c()
  
  for(i in 1:ncol(df)) {
    
    if(sum(is.na(df[, i])) / nrow(df) > 0.4) {
      
      cols_pos <- c(cols_pos, i)
      
    } else {
      
      next
      
    }
  }
  return(cols_pos)
}

cols_pos <- remove_noisy_vars(mydata)

mydata <- mydata[, -cols_pos]

noisy_c <- get_noisy_countries(mydata)

noisy_c
##  [1] "ABW" "AIA" "AND" "ANT" "ASM" "BES" "BMU" "CHI" "CIV" "COK" "CUW"
## [12] "CYM" "DMA" "ESH" "FLK" "FRO" "GIB" "GLP" "GRL" "GUF" "GUM" "IMN"
## [23] "KNA" "LIE" "MAC" "MAF" "MCO" "MHL" "MNP" "MSR" "MTQ" "MYT" "NCL"
## [34] "NIU" "NRU" "PLW" "PRI" "PRK" "PYF" "REU" "SHN" "SMR" "SPM" "SXM"
## [45] "TCA" "TKL" "TUV" "TWN" "VAT" "VGB" "VIR" "WLF" "XKX"
library(VIM)

aggr_plot <- aggr(mydata, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(mydata), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))

## 
##  Variables sorted by number of missings: 
##                    Variable      Count
##                     acc_own 0.39495798
##                   education 0.36974790
##                   inflation 0.31092437
##                    refugees 0.31092437
##                matern_leave 0.25210084
##                         BoP 0.23949580
##  deaths_household_pollution 0.23529412
##                 gni_per_cap 0.22689076
##                productivity 0.22268908
##               lab_force_fem 0.21848739
##       fem_to_male_part_rate 0.21428571
##           contr_fam_workers 0.21428571
##                 gdp_per_cap 0.21008403
##            clean_fuels_2016 0.20588235
##                 vulner_empl 0.20588235
##                  credit_inf 0.19747899
##                      health 0.19327731
##         sanitation_services 0.18487395
##                       water 0.18487395
##                     emp_agr 0.18067227
##               inf_mort_rate 0.17647059
##                    sci_jour 0.16806723
##                     eco_agr 0.12605042
##                         gdp 0.11764706
##                     eco_ind 0.11764706
##                 cell_subscr 0.11344538
##                         co2 0.11344538
##                    life_exp 0.10924370
##                  int_tr_bal 0.10924370
##                   fert_rate 0.10504202
##            electricity_2016 0.10084034
##              internet_usage 0.03361345
##                     urban_p 0.03361345
##            urban_pop_growth 0.03361345
##                    renew_en 0.02521008
##                renew_electr 0.02521008
##                country_code 0.00000000
mydata_df <- mydata

#mydata_numeric <-  data.frame(sapply(mydata[ , -c(1,2)], as.numeric))

mydata_df <- mydata_df[!duplicated(mydata_df$country_code), ]

mydata_df[ , 2:ncol(mydata_df)] <- sapply(mydata_df[ , 2:ncol(mydata_df)], as.numeric)

Ahora vamos a ver cómo quedaría el missing data plot usando sólo aquellos países con menos de 30% de valores ausentes.

require(lava)

mydata2 <- mydata_df[mydata_df$country_code %ni% noisy_c, ]

aggr_plot <- aggr(mydata2, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(mydata), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))

## 
##  Variables sorted by number of missings: 
##                    Variable       Count
##                   education 0.256830601
##                     acc_own 0.229508197
##                    refugees 0.153005464
##                   inflation 0.147540984
##                         BoP 0.092896175
##                matern_leave 0.065573770
##                productivity 0.049180328
##                 gni_per_cap 0.043715847
##                 gdp_per_cap 0.038251366
##               lab_force_fem 0.038251366
##       fem_to_male_part_rate 0.032786885
##           contr_fam_workers 0.032786885
##            clean_fuels_2016 0.027322404
##  deaths_household_pollution 0.027322404
##                     emp_agr 0.027322404
##                 vulner_empl 0.027322404
##                      health 0.021857923
##         sanitation_services 0.016393443
##                       water 0.016393443
##                         co2 0.010928962
##                    sci_jour 0.010928962
##                  credit_inf 0.010928962
##               inf_mort_rate 0.010928962
##            urban_pop_growth 0.010928962
##                     eco_agr 0.005464481
##                country_code 0.000000000
##              internet_usage 0.000000000
##                 cell_subscr 0.000000000
##                         gdp 0.000000000
##                     urban_p 0.000000000
##                    renew_en 0.000000000
##                renew_electr 0.000000000
##            electricity_2016 0.000000000
##                   fert_rate 0.000000000
##                    life_exp 0.000000000
##                  int_tr_bal 0.000000000
##                     eco_ind 0.000000000

Eliminamos algunas variables que hemos visto que, o no contienen información relevante para el análisis (hemos visto al hacer la matriz de correlaciones que no tienen demasiada correlación con el resto de variables) o hemos tenido que descartar por contener demasiados valores ausentes.

mydata2$internet_usage <- NULL

mydata2$matern_leave <- NULL

mydata2$gdp <- NULL

mydata2$BoP <- NULL

#mydata2$cell_subscr <- NULL

mydata2$contr_fam_workers <- NULL

mydata2$refugees <- NULL

#mydata2$health <- NULL

mydata2$education <- NULL

#mydata2$co2 <- NULL

#mydata2$renew_electr <- NULL

mydata2$int_forgiven <- NULL

mydata2$imf_charges <- NULL

mydata2$int_prop <- NULL

mydata2$food_deficit <- NULL

mydata2$patent_nr <- NULL

mydata2$patent_r <- NULL

#mydata2$sci_jour <- NULL

mydata2$trademark.app <- NULL

mydata2$hightech_exp <- NULL

mydata2$int_tr_bal <- NULL

mydata2$eco_ind <- NULL

mydata2$fem_to_male_part_rate <- NULL

mydata2$lab_force_fem <- NULL

mydata2$co2 <- NULL

mydata2$health <- NULL

mydata2$sci_jour <- NULL

mydata2$cell_subscr <- NULL

Vemos que mejora considerablemente, por lo que probablemente merezca la pena quedarnos solo con 180 observaciones.

if(!require(mice)) {
  
  install.packages("mice")
  
}
## Loading required package: mice
## Loading required package: lattice
## 
## Attaching package: 'mice'
## The following object is masked from 'package:tidyr':
## 
##     complete
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
countries_and_regions <- country_profile_vars[ , c("country_code","country", "Region")]

names(countries_and_regions) <- c("country_code", "country", "region")

mydata2 <- merge(mydata2, countries_and_regions, on = "country_code", all.x = T)

mydata2$region <- as.factor(mydata2$region)

mydata2$country_code <- NULL

countries <- mydata2$country

mydata2$country <- NULL

mydata_imputed <- mice(mydata2, m=1, method='cart', printFlag=FALSE)
## Warning: Number of logged events: 35
mydata_complete <- mice::complete(mydata_imputed)

sum(sapply(mydata_complete, function(x) { sum(is.na(x)) }))
## [1] 0

Las variables con las que nos quedamos son las que vienen a continuación. No explicamos todas porque la mayoría son autoexplicativas con el nombre, y por no ser demasiado redundantes, pues el significado de las variables se comenta de forma tangencial a lo largo de todo el análisis.

  1. Gdp Per Capita: la sacamos de WB
  2. Inflation: la sacamos de WB. Medida como índice de precios al consumo.
  3. Gni Per Capita: la sacamos de WB
  4. Urban Population: porcentaje de población urbana. La sacamos de UN.
  5. Employment Agriculture: porcentaje de empleo en agricultura. UN
  6. Economy Agriculture: Importancia de la agricultura en la economía. UN.
  7. Urban Population Growth: de 2016 a 2017. UN.
  8. Credit Information: Del 1 al 8 la información de crédito que tiene la gente de un país. WB.
  9. Productivity: GDP/Hour worked. WB.
  10. Renewable Energy: Uso de energía renovable. WB.
  11. Renewable Electricity: Porcentaje de la energía eléctrica que es renovable. WB.
  12. Electricity 2016: WB.
  13. Clean Fuels 2016: Acceso a “clean fuels”. WB.
  14. Infant Mortality Rate: WHO.
  15. Sanitation Services: WHO.
  16. Deaths for Household Pollution: WHO.
  17. Water: WHO.
  18. Fertility Rate: UN.
  19. Account Ownership: WB.
  20. Life Expectancy: UN.
  21. Vulnerable Employment: WB.

2. Análisis Descriptivo

2.1. Individual

if(!require(fitdistrplus)){
  
  install.packages("fitdistrplus")
  
}
## Loading required package: fitdistrplus
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: survival
## Loading required package: npsurv
## Loading required package: lsei
regions <- mydata_complete$region

mydata_complete$region <- NULL #ya no lo necesitamos

datos_brutos <- data.frame(apply(mydata_complete, 2, scale))

  
for(i in 1:ncol(datos_brutos)) { 
  
  require(ggplot2)
  
  if(is.numeric(datos_brutos[ , i])) {
    
    hist(datos_brutos[ , i], main = names(datos_brutos)[i], col = "blue")
    
    boxplot(datos_brutos[ , i], col = "green", main = names(datos_brutos)[i])
    
    if(is.integer(datos_brutos[ , i])) {
      
      descdist(datos_brutos[ , i], discrete = T) #si todos los elementos del vector son enteros, es una variable discreta que toma suficientes valores. 
      
    } else {
      
      descdist(datos_brutos[ , i], discrete = F) #sino, es una variable continua
      
    }
    
  } else if(is.factor(datos_brutos[ , i])) {
    
    counts <- table(datos_brutos[ , i])
    
    barplot(counts, main = names(datos_brutos)[i])
    
    descdist(as.numeric(datos_brutos[ , i]), discrete = T)
    

  } else {
    
    next
    
  }
  
}
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:lava':
## 
##     vars

Para no alargar demasiado unas explicaciones que tampoco tienen una relevancia crucial en todos los casos en la consecución del análisis multivariante, en esta sección nos centraremos principalmente en visualizar la distribución de todas las variables (escaladas), de tal forma que nos podremos hacer una idea de los outliers que tienen (podrían sesgar el análisis posterior). Para ello, dibujaremos el boxplot e histograma de cada una de las variables, así como, ayudándonos de la librería fitdistrplus, la distribución a la que más se asemejan las mismas. Esto nos permite, de forma sencilla y visual, la curtosis y asimetría que tienen cada una de nuestras variables, y poder compararlas respecto a las esperables para cada una de las distribuciones que incluye (ver leyenda de los gráficos). El gdp per capita se distribuye más o menos como una beta; podemos ver que tiene una fuerte asimetría a la derecha, y, viendo el boxplot, hay varias observaciones que serían outliers en una normal. La variable inflación presenta rarezas, ya que casi todos los países están en torno al 0, mientras que vemos incluso una observación en 12.31 (1231% de inflación). Esto muestra que la mayoría de países tienen una inflación tremendamente similar, siendo muy pequeño el IQR, mientras que algunos países tienen especiales problemas de inflación y contienen outliers. Estas observaciones pueden resultar problemáticas a la hora de aplicar una técnica de reducción de la dimensión como Principal Component Analysis (PCA) y Factor Analysis (FA). Debido a la enorme magnitud de esa observación, vamos a eliminarla. Lo comentado para el gdp per capita parece cumplirse también para otras variables como el gni per capita, productivity, eco_agr (esta variable representa la importancia de la agricultura como sector económico dentro de la economía doméstica). Hay otras variables que también parecen comportarse como una beta, aunque con distintos grados tanto de curtosis como de asimetría. Algo que observamos en muchas variables es que parece haber diferentes “montañas” dentro de la distribución, mostrando que puede haber distintos grupos dentro de los países, dentro de los cuales los datos sí fluctúan más o menos al rededor de una media pero que están separadas por “valles”, es decir países intermedios, siendo el número de estos bastante escaso. Un ejemplo de eso es la variable credit_inf, es decir información sobre crédito de la población. Podemos ver que hay dos medias, una para el mundo desarrollado y otra para los países en vías de desarrollo. Hay muchos países en el que el nivel de conocimiento de servicios crediticios etc es muy inferior al que podemos tener en determinados países de Europa o en América del Norte. Muy relacionada con esta variable está acc_own, es decir el porcentaje de la población que tiene cuenta bancaria o similar (explicación más profunda más arriba). Y también en esta variable podemos ver estas dos “montañas” que antes comentaba, esa separación de la muestra en dos grupos claramente diferenciados. Fert_rate es otra variable que tiene un comportamiento de este tipo, aunque mucho menos pronunciado. Las variables que se comportan de esta forma son, a priori, las que mejor nos van a ayudar a distinguir el nivel de spread, o el posicionamiento con respecto al evento spread (es decir qué tal posicionado está el país con respecto a los demás en términos de calidad de vida derivada del progreso económico y tecnológico). El resto de variables siguen distribuciones o bien similares a la uniforme, o bien similares a la beta con distintos grados de curtosis y también de asimetría (además, tenemos asimetrías positivas y negativas dentro de la muestra, no todas las variables la tienen en la misma dirección). Aparte de los ya comentados outliers de inflation, no observamos ningún otro desajuste excesivo en los datos, por lo que podemos proceder a continuar con la parte más importante del análisis descriptivo, que es la multivariante. Más que cómo están distribuidas nuestras variables, nos interesa saber cómo se distribuyen unas con respecto a las otras; queremos ver si se adivinan correlaciones lineales entre las variables. En definitiva el objetivo es tener un criterio previo que nos defina qué variables pueden o pueden no ser adecuadas y significativas para el análisis posterior, con el fin de eliminar ruido de la muestra y tener una idea previa de la cantidad de información común entre las variables. Cuanta más relación contengan, sin llegar a tener información redundante (necesitamos algo de multicolinearidad para que se puedan extraer factores significativos, pero ésta no puede ser perfecta porque no podríamos utilizar el método necesario para la obtención de factores – matriz no invertible).

datos_brutos <- datos_brutos[!datos_brutos$inflation == max(datos_brutos$inflation), ]

2.2. Multivariante

require(GGally)

#datos_brutos <- sapply(mydata[ , -1], as.numeric)

#datos_brutos$inflation <- NULL

#datos_brutos$BoP <- NULL

cormatrix <- cor(datos_brutos)

ggcorr(cormatrix, label = T) 

ggpairs(datos_brutos)

En los gráficos superiores se pueden apreciar la matriz de correlación de nuestras variables. En primer lugar, tenemos este gráfico de la matriz de correlación simple, en el que, como podemos ver, las variables tienen una correlación tremendamente alta unas con otras (tanto positiva como negativa). Aunque en muchos casos pone 1, es porque lo aproxima, y en el siguiente gráfico podemos ver que no llegan nunca a ser 1. Nos hemos asegurado antes de continuar de que no teníamos ningún caso de multicolinearidad perfecta, por lo comentado en el apartado anterior. En el segundo de estos gráficos, además, tenemos el scatterplot de unas variables con respecto a otras, y en la diagonal la pdf, en la que podemos ver cómo se distribuye cada variable, similar a los histogramas realizados antes. Podemos ver que muchas variables tienen relación lineal con las demás.

boxplot(datos_brutos, col="blue", las=2)

library(igraph)
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:lava':
## 
##     compare, path
## The following object is masked from 'package:tidyr':
## 
##     crossing
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
R <- cor(datos_brutos)

# visualize size of correlations
heatmap(R, symm = TRUE)

# better to reorder colors
heatmap(sqrt(2*(1-R)), symm = TRUE)  # using a distance through correlation

# Correlations induce some network structure 
graph.f<-graph.adjacency(R, weighted=TRUE, mode="undirected", diag=F)
plot.igraph(graph.f, edge.width=2,vertex.size=20)

# remove correlations lower than than percentile 20%
g=delete.edges(graph.f, E(graph.f)[ abs(weight) < quantile(abs(R), probs=0.2)])
plot.igraph(g, edge.width=2,vertex.size=20)

# But the relevant multivatiate information is in the inverse of R
iR <- solve(R)
igraph.f<-graph.adjacency(iR, weighted=TRUE, mode="lower", diag=F)
# remove partial correlations lower than perentile 20%
ig <- delete.edges(igraph.f, E(igraph.f)[ abs(weight) < quantile(abs(iR), probs=0.2) ])
plot.igraph(ig, edge.width=2,vertex.size=20)

En este chunk de arriba vemos diferentes formas de visualizar la información que contiene nuestro set de datos, en concreto la información relativa a cómo se relacionan unas variables con otras.

2.3. Inferencia

En este apartado realizaremos algunos test de hipótesis; queremos ver si nuestros datos siguen una distribución normal, entre otras cosas.

datos_brutos = as.matrix(datos_brutos)
n = dim(datos_brutos)[1]
p = dim(datos_brutos)[2]
n
## [1] 182
p
## [1] 21
# Check whether mu=0 (vector hypothesis)
mu0 = rep(0,1,p)
Sigma0 <- (1/n) * t(datos_brutos) %*% datos_brutos - mu0 %*% t(mu0)
Sigma0
##                            gdp_per_cap   inflation gni_per_cap     urban_p
## gdp_per_cap                 0.99743869 -0.10779893  0.98791782  0.59331094
## inflation                  -0.10779893  0.14310958 -0.10909547 -0.08746353
## gni_per_cap                 0.98791782 -0.10909547  0.99793090  0.60364108
## urban_p                     0.59331094 -0.08746353  0.60364108  0.98543922
## productivity                0.88209569 -0.11348328  0.87177861  0.66279367
## credit_inf                  0.25899746 -0.05042613  0.28207298  0.41305309
## renew_en                   -0.31505805  0.11550050 -0.31936175 -0.48033537
## renew_electr               -0.02461287  0.09524137 -0.02939236 -0.17337863
## electricity_2016            0.38913598 -0.12384488  0.39711980  0.53094788
## clean_fuels_2016            0.52316039 -0.09242545  0.53158870  0.64111485
## inf_mort_rate              -0.52237889  0.13830966 -0.53320946 -0.55622284
## sanitation_services         0.50727895 -0.10420559  0.51792597  0.57458375
## deaths_household_pollution -0.54700681  0.08785447 -0.55532721 -0.61253540
## water                       0.46378713 -0.13528488  0.47426166  0.59812316
## fert_rate                  -0.47746874  0.12291346 -0.48791983 -0.52964720
## life_exp                    0.63628922 -0.14348935  0.64967537  0.62261185
## urban_pop_growth           -0.30348990  0.07348143 -0.31300274 -0.30968505
## emp_agr                    -0.57247374  0.10621856 -0.58193154 -0.69701655
## eco_agr                    -0.51314458  0.10637209 -0.52011440 -0.59107134
## vulner_empl                -0.60612344  0.10856046 -0.61768529 -0.64416840
## acc_own                     0.70775770 -0.12814266  0.71954776  0.52866905
##                            productivity  credit_inf   renew_en
## gdp_per_cap                   0.8820957  0.25899746 -0.3150581
## inflation                    -0.1134833 -0.05042613  0.1155005
## gni_per_cap                   0.8717786  0.28207298 -0.3193618
## urban_p                       0.6627937  0.41305309 -0.4803354
## productivity                  0.9947228  0.33485030 -0.5035063
## credit_inf                    0.3348503  0.98731418 -0.3088054
## renew_en                     -0.5035063 -0.30880539  0.9996920
## renew_electr                 -0.1844359  0.06198460  0.5586055
## electricity_2016              0.5193091  0.48992052 -0.7219390
## clean_fuels_2016              0.6426168  0.49230892 -0.7289692
## inf_mort_rate                -0.6015568 -0.53903195  0.6212961
## sanitation_services           0.6225482  0.51325062 -0.7252856
## deaths_household_pollution   -0.6312588 -0.46800456  0.6701490
## water                         0.5673666  0.47216510 -0.6821997
## fert_rate                    -0.5604317 -0.54391887  0.6201322
## life_exp                      0.6750055  0.49890395 -0.6317957
## urban_pop_growth             -0.3571955 -0.32136303  0.5096508
## emp_agr                      -0.6966648 -0.40254688  0.7205614
## eco_agr                      -0.6337788 -0.45775220  0.6663125
## vulner_empl                  -0.7287212 -0.43582621  0.6937656
## acc_own                       0.7259118  0.45017963 -0.4540027
##                            renew_electr electricity_2016 clean_fuels_2016
## gdp_per_cap                 -0.02461287        0.3891360       0.52316039
## inflation                    0.09524137       -0.1238449      -0.09242545
## gni_per_cap                 -0.02939236        0.3971198       0.53158870
## urban_p                     -0.17337863        0.5309479       0.64111485
## productivity                -0.18443594        0.5193091       0.64261682
## credit_inf                   0.06198460        0.4899205       0.49230892
## renew_en                     0.55860548       -0.7219390      -0.72896921
## renew_electr                 0.99446396       -0.2009990      -0.21720152
## electricity_2016            -0.20099898        0.9638105       0.83557168
## clean_fuels_2016            -0.21720152        0.8355717       0.98421657
## inf_mort_rate                0.18683577       -0.7766352      -0.78362738
## sanitation_services         -0.19072477        0.8398051       0.86584425
## deaths_household_pollution   0.18878617       -0.7020436      -0.87447132
## water                       -0.22503654        0.8393339       0.78350552
## fert_rate                    0.14273061       -0.7847359      -0.76916617
## life_exp                    -0.17072089        0.7870422       0.77999648
## urban_pop_growth             0.14962094       -0.5945882      -0.61503664
## emp_agr                      0.32390290       -0.7422659      -0.82466812
## eco_agr                      0.19220274       -0.6784794      -0.74405649
## vulner_empl                  0.26598501       -0.7307746      -0.84512504
## acc_own                     -0.13166022        0.5272506       0.60586083
##                            inf_mort_rate sanitation_services
## gdp_per_cap                   -0.5223789           0.5072789
## inflation                      0.1383097          -0.1042056
## gni_per_cap                   -0.5332095           0.5179260
## urban_p                       -0.5562228           0.5745838
## productivity                  -0.6015568           0.6225482
## credit_inf                    -0.5390319           0.5132506
## renew_en                       0.6212961          -0.7252856
## renew_electr                   0.1868358          -0.1907248
## electricity_2016              -0.7766352           0.8398051
## clean_fuels_2016              -0.7836274           0.8658443
## inf_mort_rate                  0.9777406          -0.8253493
## sanitation_services           -0.8253493           0.9749646
## deaths_household_pollution     0.7803206          -0.7905230
## water                         -0.7905601           0.8450177
## fert_rate                      0.8309591          -0.8153358
## life_exp                      -0.9041731           0.8323763
## urban_pop_growth               0.5805201          -0.6309816
## emp_agr                        0.7374044          -0.8024338
## eco_agr                        0.6977185          -0.7190240
## vulner_empl                    0.7712448          -0.8219014
## acc_own                       -0.6912545           0.6195550
##                            deaths_household_pollution      water
## gdp_per_cap                               -0.54700681  0.4637871
## inflation                                  0.08785447 -0.1352849
## gni_per_cap                               -0.55532721  0.4742617
## urban_p                                   -0.61253540  0.5981232
## productivity                              -0.63125881  0.5673666
## credit_inf                                -0.46800456  0.4721651
## renew_en                                   0.67014901 -0.6821997
## renew_electr                               0.18878617 -0.2250365
## electricity_2016                          -0.70204365  0.8393339
## clean_fuels_2016                          -0.87447132  0.7835055
## inf_mort_rate                              0.78032055 -0.7905601
## sanitation_services                       -0.79052297  0.8450177
## deaths_household_pollution                 0.98807701 -0.6957409
## water                                     -0.69574090  0.9779317
## fert_rate                                  0.67233025 -0.8176361
## life_exp                                  -0.77115633  0.7972351
## urban_pop_growth                           0.48726028 -0.6298389
## emp_agr                                    0.76058708 -0.7688848
## eco_agr                                    0.70599806 -0.6559515
## vulner_empl                                0.79936960 -0.7424831
## acc_own                                   -0.60697084  0.6063777
##                             fert_rate   life_exp urban_pop_growth
## gdp_per_cap                -0.4774687  0.6362892      -0.30348990
## inflation                   0.1229135 -0.1434894       0.07348143
## gni_per_cap                -0.4879198  0.6496754      -0.31300274
## urban_p                    -0.5296472  0.6226118      -0.30968505
## productivity               -0.5604317  0.6750055      -0.35719553
## credit_inf                 -0.5439189  0.4989040      -0.32136303
## renew_en                    0.6201322 -0.6317957       0.50965075
## renew_electr                0.1427306 -0.1707209       0.14962094
## electricity_2016           -0.7847359  0.7870422      -0.59458823
## clean_fuels_2016           -0.7691662  0.7799965      -0.61503664
## inf_mort_rate               0.8309591 -0.9041731       0.58052005
## sanitation_services        -0.8153358  0.8323763      -0.63098164
## deaths_household_pollution  0.6723303 -0.7711563       0.48726028
## water                      -0.8176361  0.7972351      -0.62983885
## fert_rate                   0.9860985 -0.8225073       0.69205254
## life_exp                   -0.8225073  0.9792918      -0.55331598
## urban_pop_growth            0.6920525 -0.5533160       0.98308739
## emp_agr                     0.7156742 -0.7493805       0.60936854
## eco_agr                     0.6637815 -0.6713800       0.50745172
## vulner_empl                 0.7432691 -0.7550921       0.62703060
## acc_own                    -0.6857094  0.6957487      -0.49989456
##                               emp_agr    eco_agr vulner_empl    acc_own
## gdp_per_cap                -0.5724737 -0.5131446  -0.6061234  0.7077577
## inflation                   0.1062186  0.1063721   0.1085605 -0.1281427
## gni_per_cap                -0.5819315 -0.5201144  -0.6176853  0.7195478
## urban_p                    -0.6970165 -0.5910713  -0.6441684  0.5286691
## productivity               -0.6966648 -0.6337788  -0.7287212  0.7259118
## credit_inf                 -0.4025469 -0.4577522  -0.4358262  0.4501796
## renew_en                    0.7205614  0.6663125   0.6937656 -0.4540027
## renew_electr                0.3239029  0.1922027   0.2659850 -0.1316602
## electricity_2016           -0.7422659 -0.6784794  -0.7307746  0.5272506
## clean_fuels_2016           -0.8246681 -0.7440565  -0.8451250  0.6058608
## inf_mort_rate               0.7374044  0.6977185   0.7712448 -0.6912545
## sanitation_services        -0.8024338 -0.7190240  -0.8219014  0.6195550
## deaths_household_pollution  0.7605871  0.7059981   0.7993696 -0.6069708
## water                      -0.7688848 -0.6559515  -0.7424831  0.6063777
## fert_rate                   0.7156742  0.6637815   0.7432691 -0.6857094
## life_exp                   -0.7493805 -0.6713800  -0.7550921  0.6957487
## urban_pop_growth            0.6093685  0.5074517   0.6270306 -0.4998946
## emp_agr                     0.9765394  0.7921508   0.8885272 -0.6508890
## eco_agr                     0.7921508  0.9978747   0.7816005 -0.6176567
## vulner_empl                 0.8885272  0.7816005   0.9986436 -0.7075213
## acc_own                    -0.6508890 -0.6176567  -0.7075213  0.9793095
S = cov(datos_brutos)  
S
##                            gdp_per_cap   inflation gni_per_cap     urban_p
## gdp_per_cap                 1.00293525 -0.10813568  0.99336321  0.59655516
## inflation                  -0.10813568  0.13916604 -0.10946558 -0.08732963
## gni_per_cap                 0.99336321 -0.10946558  1.00343290  0.60694578
## urban_p                     0.59655516 -0.08732963  0.60694578  0.99080319
## productivity                0.88694883 -0.11373873  0.87657681  0.66640708
## credit_inf                  0.26039690 -0.05012870  0.28360309  0.41526006
## renew_en                   -0.31679380  0.11604887 -0.32112177 -0.48297745
## renew_electr               -0.02476966  0.09614809 -0.02957345 -0.17438613
## electricity_2016            0.39123271 -0.12355619  0.39926602  0.53375447
## clean_fuels_2016            0.52601565 -0.09229357  0.53449408  0.64457316
## inf_mort_rate              -0.52522325  0.13831077 -0.53611787 -0.55919643
## sanitation_services         0.51003735 -0.10397210  0.52074767  0.57765276
## deaths_household_pollution -0.54999842  0.08778141 -0.55836787 -0.61584678
## water                       0.46630796 -0.13527256  0.47684455  0.60132867
## fert_rate                  -0.48007372  0.12298955 -0.49058589 -0.53249483
## life_exp                    0.63976439 -0.14354615  0.65322858  0.62595576
## urban_pop_growth           -0.30513028  0.07322231 -0.31469936 -0.31130931
## emp_agr                    -0.57559375  0.10602206 -0.58510814 -0.70076535
## eco_agr                    -0.51599253  0.10719555 -0.52299955 -0.59436766
## vulner_empl                -0.60946189  0.10897189 -0.62108866 -0.64770279
## acc_own                     0.71162774 -0.12811498  0.72348701  0.53149398
##                            productivity  credit_inf   renew_en
## gdp_per_cap                   0.8869488  0.26039690 -0.3167938
## inflation                    -0.1137387 -0.05012870  0.1160489
## gni_per_cap                   0.8765768  0.28360309 -0.3211218
## urban_p                       0.6664071  0.41526006 -0.4829775
## productivity                  1.0001894  0.33665510 -0.5062811
## credit_inf                    0.3366551  0.99269887 -0.3105006
## renew_en                     -0.5062811 -0.31050058  1.0052134
## renew_electr                 -0.1854848  0.06228075  0.5616989
## electricity_2016              0.5221019  0.49250889 -0.7259092
## clean_fuels_2016              0.6461168  0.49495068 -0.7329845
## inf_mort_rate                -0.6048204 -0.54191718  0.6247142
## sanitation_services           0.6259242  0.51598780 -0.7292774
## deaths_household_pollution   -0.6347026 -0.47052228  0.6738409
## water                         0.5704416  0.47468130 -0.6859544
## fert_rate                    -0.5634807 -0.54685058  0.6235469
## life_exp                      0.6786770  0.50157078 -0.6352723
## urban_pop_growth             -0.3591168 -0.32305759  0.5124539
## emp_agr                      -0.7004523 -0.40467558  0.7245275
## eco_agr                      -0.6372989 -0.46030991  0.6699983
## vulner_empl                  -0.7327325 -0.43821117  0.6975950
## acc_own                       0.7298646  0.45257730 -0.4564970
##                            renew_electr electricity_2016 clean_fuels_2016
## gdp_per_cap                 -0.02476966        0.3912327       0.52601565
## inflation                    0.09614809       -0.1235562      -0.09229357
## gni_per_cap                 -0.02957345        0.3992660       0.53449408
## urban_p                     -0.17438613        0.5337545       0.64457316
## productivity                -0.18548479        0.5221019       0.64611677
## credit_inf                   0.06228075        0.4925089       0.49495068
## renew_en                     0.56169892       -0.7259092      -0.73298449
## renew_electr                 0.99992765       -0.2021877      -0.21845317
## electricity_2016            -0.20218767        0.9689354       0.84005605
## clean_fuels_2016            -0.21845317        0.8400561       0.98956703
## inf_mort_rate                0.18792934       -0.7807692      -0.78785326
## sanitation_services         -0.19184354        0.8442786       0.87051810
## deaths_household_pollution   0.18987407       -0.7058076      -0.87922686
## water                       -0.22634090        0.8438150       0.78773117
## fert_rate                    0.14356765       -0.7889475      -0.77333387
## life_exp                    -0.17172326        0.7912392       0.78420597
## urban_pop_growth             0.15050103       -0.5977366      -0.61834437
## emp_agr                      0.32575539       -0.7462058      -0.82911799
## eco_agr                      0.19324569       -0.6822764      -0.74819930
## vulner_empl                  0.26746968       -0.7347733      -0.84976867
## acc_own                     -0.13244676        0.5300124       0.60910828
##                            inf_mort_rate sanitation_services
## gdp_per_cap                   -0.5252232           0.5100374
## inflation                      0.1383108          -0.1039721
## gni_per_cap                   -0.5361179           0.5207477
## urban_p                       -0.5591964           0.5776528
## productivity                  -0.6048204           0.6259242
## credit_inf                    -0.5419172           0.5159878
## renew_en                       0.6247142          -0.7292774
## renew_electr                   0.1879293          -0.1918435
## electricity_2016              -0.7807692           0.8442786
## clean_fuels_2016              -0.7878533           0.8705181
## inf_mort_rate                  0.9830195          -0.8297788
## sanitation_services           -0.8297788           0.9802128
## deaths_household_pollution     0.7845417          -0.7947950
## water                         -0.7948054           0.8495564
## fert_rate                      0.8354528          -0.8197373
## life_exp                      -0.9090499           0.8368492
## urban_pop_growth               0.5836201          -0.6343540
## emp_agr                        0.7413522          -0.8067332
## eco_agr                        0.7016113          -0.7230368
## vulner_empl                    0.7754755          -0.8264101
## acc_own                       -0.6949551           0.6228522
##                            deaths_household_pollution      water
## gdp_per_cap                               -0.54999842  0.4663080
## inflation                                  0.08778141 -0.1352726
## gni_per_cap                               -0.55836787  0.4768446
## urban_p                                   -0.61584678  0.6013287
## productivity                              -0.63470260  0.5704416
## credit_inf                                -0.47052228  0.4746813
## renew_en                                   0.67384090 -0.6859544
## renew_electr                               0.18987407 -0.2263409
## electricity_2016                          -0.70580758  0.8438150
## clean_fuels_2016                          -0.87922686  0.7877312
## inf_mort_rate                              0.78454171 -0.7948054
## sanitation_services                       -0.79479505  0.8495564
## deaths_household_pollution                 0.99347013 -0.6994952
## water                                     -0.69949516  0.9832127
## fert_rate                                  0.67597366 -0.8220567
## life_exp                                  -0.77533005  0.8015216
## urban_pop_growth                           0.48987387 -0.6332119
## emp_agr                                    0.76469682 -0.7730071
## eco_agr                                    0.70992641 -0.6596134
## vulner_empl                                0.80376379 -0.7465550
## acc_own                                   -0.61023749  0.6096098
##                             fert_rate   life_exp urban_pop_growth
## gdp_per_cap                -0.4800737  0.6397644      -0.30513028
## inflation                   0.1229895 -0.1435461       0.07322231
## gni_per_cap                -0.4905859  0.6532286      -0.31469936
## urban_p                    -0.5324948  0.6259558      -0.31130931
## productivity               -0.5634807  0.6786770      -0.35911679
## credit_inf                 -0.5468506  0.5015708      -0.32305759
## renew_en                    0.6235469 -0.6352723       0.51245389
## renew_electr                0.1435676 -0.1717233       0.15050103
## electricity_2016           -0.7889475  0.7912392      -0.59773657
## clean_fuels_2016           -0.7733339  0.7842060      -0.61834437
## inf_mort_rate               0.8354528 -0.9090499       0.58362015
## sanitation_services        -0.8197373  0.8368492      -0.63435404
## deaths_household_pollution  0.6759737 -0.7753300       0.48987387
## water                      -0.8220567  0.8015216      -0.63321189
## fert_rate                   0.9914698 -0.8269578       0.69579132
## life_exp                   -0.8269578  0.9845879      -0.55626958
## urban_pop_growth            0.6957913 -0.5562696       0.98842537
## emp_agr                     0.7195284 -0.7533990       0.61262517
## eco_agr                     0.6674789 -0.6751259       0.51028844
## vulner_empl                 0.7473516 -0.7592346       0.63046839
## acc_own                    -0.6894042  0.6994783      -0.50255305
##                               emp_agr    eco_agr vulner_empl    acc_own
## gdp_per_cap                -0.5755938 -0.5159925  -0.6094619  0.7116277
## inflation                   0.1060221  0.1071956   0.1089719 -0.1281150
## gni_per_cap                -0.5851081 -0.5229995  -0.6210887  0.7234870
## urban_p                    -0.7007654 -0.5943677  -0.6477028  0.5314940
## productivity               -0.7004523 -0.6372989  -0.7327325  0.7298646
## credit_inf                 -0.4046756 -0.4603099  -0.4382112  0.4525773
## renew_en                    0.7245275  0.6699983   0.6975950 -0.4564970
## renew_electr                0.3257554  0.1932457   0.2674697 -0.1324468
## electricity_2016           -0.7462058 -0.6822764  -0.7347733  0.5300124
## clean_fuels_2016           -0.8291180 -0.7481993  -0.8497687  0.6091083
## inf_mort_rate               0.7413522  0.7016113   0.7754755 -0.6949551
## sanitation_services        -0.8067332 -0.7230368  -0.8264101  0.6228522
## deaths_household_pollution  0.7646968  0.7099264   0.8037638 -0.6102375
## water                      -0.7730071 -0.6596134  -0.7465550  0.6096098
## fert_rate                   0.7195284  0.6674789   0.7473516 -0.6894042
## life_exp                   -0.7533990 -0.6751259  -0.7592346  0.6994783
## urban_pop_growth            0.6126252  0.5102884   0.6304684 -0.5025531
## emp_agr                     0.9818050  0.7965663   0.8934050 -0.6543633
## eco_agr                     0.7965663  1.0033760   0.7859281 -0.6211058
## vulner_empl                 0.8934050  0.7859281   1.0041535 -0.7114010
## acc_own                    -0.6543633 -0.6211058  -0.7114010  0.9846057
lambda <- n * log(det(Sigma0)/det(S))
lambda
## [1] -12.50472
pvalue <- 1 - pchisq(lambda,p)
pvalue
## [1] 1
#  Accept mu=0


# Check whether Sigma=diag (independence under normality)
R = cor(datos_brutos)
lambda <- -n * log(det(R))
lambda
## [1] 5248.799
pvalue <- 1 - pchisq(lambda,p*(p-1)/2)
pvalue
## [1] 0

Como podemos ver arriba, después de hacer inferencia podemos decir que aceptamos la H0 de que \(\mu = 0\), mientras que descartamos la hipótesis de que nuestras variables estén incorreladas.

S <- cov(datos_brutos)
muhat = colMeans(datos_brutos)
dM2 <- mahalanobis(datos_brutos,muhat,S)  # Mahalanobis distances (square)
dM2
##          1          2          3          4          5          6 
##  28.994962  68.555443  21.388509   8.404907   7.558137  13.180230 
##          7          8          9         10         11         12 
##  18.374634   9.696618   7.720667  21.828417  40.022056   7.125166 
##         13         14         15         16         17         18 
##  18.604817  17.877130  24.750341   6.219358   6.096611  22.091780 
##         19         20         21         22         23         24 
##  23.295159  13.033114  15.661010  14.348899   8.837935  15.919010 
##         25         26         27         28         29         30 
##  31.224956  25.371622  15.826720  32.180285   6.444384  39.388383 
##         31         32         33         34         35         36 
##   6.182840  20.168302  15.728870  51.061244  25.523341  11.542334 
##         37         38         39         40         41         42 
##  23.726673   5.609281  15.090946  15.153698   4.929121   2.919185 
##         43         44         45         46         47         48 
##   5.552209  32.563770  11.065248   9.515942  13.065805   9.020102 
##         49         50         51         52         53         54 
##  43.520519  59.652245   6.078852   7.027677  28.643175   7.920348 
##         55         56         57         58         59         60 
##  25.020633   5.560876  28.249454  42.183325   6.815805  27.789687 
##         61         62         63         64         65         66 
##  30.990331  15.348570  21.747149  23.384304  53.329767   6.557230 
##         67         68         69         70         71         72 
##  13.205011  16.021994  13.215927  10.896540   9.845490   7.955181 
##         73         74         75         76         77         78 
##  34.236913   4.782653   7.203811  22.786987  29.744670  20.300150 
##         79         80         81         82         83         84 
##  22.821519  35.986102  12.566905   8.842708   6.403739  18.261907 
##         85         86         87         88         89         90 
##   8.159629  11.226229  30.453753  13.543354  20.859019  33.518235 
##         91         92         93         94         95         96 
##  67.035594  24.947342  37.321071  14.566408  56.126314  22.700155 
##         97         98         99        100        101        102 
##  22.080959  24.545888  41.604811   6.617250 140.925812   8.102901 
##        103        104        105        106        107        108 
##  14.297767  15.066632  21.627115  21.104485   8.191407  11.423334 
##        109        110        111        112        113        114 
##  21.665455   9.577928  31.399360  19.964891  23.727328  18.253856 
##        115        116        117        118        119        120 
##  16.961657   7.941349  28.143204   7.602743  43.726340  30.849777 
##        121        122        123        124        125        126 
##  46.365370  13.530767   6.101376  30.791715  27.825926  11.319838 
##        127        128        129        130        131        132 
##  43.361454  33.042003  13.166181   9.251674  13.843996  31.744640 
##        133        134        135        136        137        138 
##   5.289577   5.845897  13.592292  20.383757  34.692212  14.075633 
##        139        140        141        142        143        144 
##  10.032621  32.351157  23.336812  21.270008  13.903470  14.511892 
##        145        146        147        148        149        151 
##  18.214657  29.745228   9.861480  46.161559  10.564650  24.408264 
##        152        153        154        155        156        157 
##  26.068018   4.663597   8.031725  12.554002  29.258330   5.160675 
##        158        159        160        161        162        163 
##  16.341235  37.561653  21.422667  14.616840  32.627211  20.300257 
##        164        165        166        167        168        169 
##  34.837545  22.561278  28.966642   6.720782  10.242489  16.788131 
##        170        171        172        173        174        175 
##  25.786257  17.064537  13.055307  13.931525  21.598127  16.581537 
##        176        177        178        179        180        181 
##  18.760829  12.186404  23.239589  42.546579  15.428137  20.548363 
##        182        183 
##  21.887457  21.878487
plot(sort(dM2,index.return=TRUE)$x,pch=19,col=grey.colors(n=n,start=0,end=.75),xlab="",ylab="",main="Mahalanobis distances (square)")

## Q-Q plot for Chi^2 data against true theoretical distribution:
qqplot(dM2, qchisq(ppoints(n), df = p),xlab="quantiles of Mahalanobis distances",
       ylab="Quantiles for Chi-square with p df",pch=19,col="darkred",
       main = expression("Q-Q plot for" ~~ {chi^2}[p]))
qqline(dM2, distribution = function(q) qchisq(q, df = p),
       prob = c(0.1, 0.6), col = "blue")

Aquí vemos que no se distribuye de forma demasiado normal. PCA a priori no necesita de ninguna suposición de normalidad para que la extracción de la varianza máxima en los menores componentes posibles (los primeros son los que mayor porcentaje de varianza explican) sea eficiente. FA utiliza la suposición de normalidad en las uniquesses del modelo, ya que \(\epsilon \sim N(0, \sigma)\). El hecho de que nuestros datos no se distribuyan de forma normal podría afectar al comportamiento de estas uniquesses, haciendo que los test que utilizamos para ver si los factores son suficientes o no dejarían de tener sentido, pues el p-value se calcula bajo la suposición de normalidad de los datos. Por este motivo en la práctica a menudo no se hace demasiado caso a ese p-value a la hora de evaluar un modelo concreto, sino que se emplea más para comparar unos modelos con otros y ver cuáles explican mejor.

3. PCA & Factor Analysis

Como tenemos 21 variables en nuestra muestra, existen \(2^{p} \approx 2^{21} = 2097152\) diferentes relaciones. El primer objetivo es, por lo tanto, ser capaces de reducir la dimensión, limpiando el ruido. En alta dimensión casi todas las observaciones son outliers; cuando aumentamos nuestros datos, especialmente el número de variables, aumenta la información pero el ruido lo hace de una forma más rápida. Por eso con cada dimensión adicional que incluyamos en nuestro set de datos estaremos incluyendo mucho más ruido que información.

Inicialmente usaremos PCA para detectar cuánta varianza común tienen nuestras variables, para después tratar de encontrar y definir las variables mencionadas al principio del estudio, en especial la variable spread. Como ya se comentó antes, el spread es en definitiva cómo se distribuye ese bounty, y es una variable latente más fácil de encontrar y de definir que el bounty, pues para estimar este necesitaríamos las variables en forma de serie temporal y, dado que para muchas variables hemos tenido que apañarnos con datos de años pasados al 2017 o con medias de los últimos años (para subsanar la ausencia de observaciones más recientes en muchos países), esto sería una tarea infactible muy probablemente. Por eso nos centraremos principalmente en encontrar el spread, es decir la creciente desigualdad causada en gran parte por el desarrollo tecnológico cuyos frutos están desigualmente distribuidos.

FA es lo mismo que PCA si le quitas las uniquesses (o si estas son 0):

PCA:

\(\Sigma = V \Lambda V^{T}\)

FA:

\(\Sigma = L L^{T} + \Psi\)

Si \(\Psi = 0\), FA es igual que PCA con \(L = V\Lambda^{1/2}\).

Para poder definir formalmente la variable spread, debemos buscarla en forma de variable latente mediante un Factor Analysis, que nos permita además ver cuánta varianza explica y deja de explicar la variable oculta spread. Sin embargo, realizamos el PCA primero para ver con cuántos factores será buena idea dotar al modelo, cuánta varianza común hay entre nuestras variables y la posibilidad de reducir la dimensión eficientemente.

3.1. Principal Component Analysis

pca <- prcomp(datos_brutos, center = T, scale = T)

screeplot(pca,main="Screeplot",col="blue",type="barplot",pch=19)

pca
## Standard deviations (1, .., p=21):
##  [1] 3.59763732 1.37042017 1.14657322 0.96253675 0.84945791 0.73966812
##  [7] 0.67354053 0.64932915 0.55332286 0.54219114 0.46175454 0.44023830
## [13] 0.37831575 0.37058149 0.35487052 0.32479451 0.30399767 0.27949216
## [19] 0.23136380 0.20584373 0.09165425
## 
## Rotation (n x k) = (21 x 21):
##                                   PC1         PC2         PC3          PC4
## gdp_per_cap                 0.1924576  0.49691736 -0.07803238  0.003498832
## inflation                  -0.1079339 -0.01727274  0.33725021  0.797897121
## gni_per_cap                 0.1952434  0.48950047 -0.06895215 -0.002237555
## urban_p                     0.2022804  0.16050288 -0.05623977  0.213630935
## productivity                0.2199064  0.35303008 -0.14848017  0.095576458
## credit_inf                  0.1532914 -0.05586961  0.42122225 -0.115073872
## renew_en                   -0.2102537  0.27506181  0.29521566 -0.142928811
## renew_electr               -0.0717081  0.25033208  0.69820706 -0.090506917
## electricity_2016            0.2393336 -0.22448831  0.08555264 -0.058065324
## clean_fuels_2016            0.2539321 -0.11149131  0.05623341  0.167900170
## inf_mort_rate              -0.2492296  0.07778861 -0.10960187  0.159272093
## sanitation_services         0.2549754 -0.13437552  0.09746871  0.034300013
## deaths_household_pollution -0.2394578  0.03309485 -0.04174147 -0.185848424
## water                       0.2439229 -0.15506465  0.05377820 -0.107947086
## fert_rate                  -0.2416519  0.12205732 -0.16985839  0.182212338
## life_exp                    0.2538097  0.01119038  0.06984229 -0.146461950
## urban_pop_growth           -0.1846428  0.20705657 -0.12273970  0.077605863
## emp_agr                    -0.2521812  0.04464787  0.10093908 -0.191426824
## eco_agr                    -0.2289405  0.03035463  0.01017343 -0.133410786
## vulner_empl                -0.2533910  0.01013596  0.04523106 -0.150156519
## acc_own                     0.2195788  0.21986228  0.01935181 -0.157304329
##                                     PC5          PC6         PC7
## gdp_per_cap                -0.114164775  0.005083875 -0.11581999
## inflation                  -0.172176012 -0.212300371 -0.31813459
## gni_per_cap                -0.100490013 -0.016017818 -0.11296041
## urban_p                     0.383880655  0.133240270  0.27499225
## productivity               -0.002476745 -0.052855788 -0.03559855
## credit_inf                  0.569439765 -0.557853795  0.02786649
## renew_en                   -0.057692902  0.008488227  0.06340036
## renew_electr               -0.101203832  0.419789420  0.26115582
## electricity_2016            0.025844979  0.240202817 -0.09373285
## clean_fuels_2016            0.017726890  0.171381240  0.00116888
## inf_mort_rate              -0.026813855 -0.029303878  0.27622622
## sanitation_services        -0.042070660  0.133868600 -0.13727783
## deaths_household_pollution -0.120129005 -0.197163678  0.12829073
## water                      -0.039932964  0.141911082 -0.09424735
## fert_rate                   0.139291134  0.105954085  0.11815650
## life_exp                    0.017329116  0.127982450 -0.34642161
## urban_pop_growth            0.621533207  0.284822827 -0.29564473
## emp_agr                     0.041977648 -0.020294054 -0.26754542
## eco_agr                    -0.092278205  0.012240581 -0.50538452
## vulner_empl                 0.095839515  0.034318499 -0.21107176
## acc_own                    -0.130596943 -0.413851902 -0.02313192
##                                     PC8         PC9         PC10
## gdp_per_cap                 0.009647005  0.21046911 -0.074492312
## inflation                   0.077825079 -0.01170776  0.213860280
## gni_per_cap                 0.019504881  0.18312985 -0.084822640
## urban_p                     0.653974934 -0.28237279 -0.005637807
## productivity               -0.026488158  0.27230759  0.036581697
## credit_inf                 -0.010439600  0.26430081 -0.216412234
## renew_en                    0.143621115 -0.36181378 -0.067124106
## renew_electr               -0.110201741  0.08718815  0.048694568
## electricity_2016            0.070694066  0.32346888  0.139484180
## clean_fuels_2016           -0.090140506  0.03667897 -0.317561885
## inf_mort_rate               0.112352458  0.37416017 -0.105206490
## sanitation_services        -0.030303241  0.15255269 -0.002555221
## deaths_household_pollution  0.336026447  0.29702805  0.486051063
## water                       0.316868790  0.12657818  0.206200369
## fert_rate                  -0.148193007  0.06991818 -0.230951167
## life_exp                    0.040733732 -0.16815145  0.100861539
## urban_pop_growth           -0.241526402 -0.05942265  0.357526995
## emp_agr                    -0.060756207  0.10328084 -0.052958568
## eco_agr                     0.405772198 -0.02144044 -0.460955589
## vulner_empl                 0.137070945  0.15146356  0.114733881
## acc_own                    -0.152099675 -0.33988491  0.237723801
##                                   PC11         PC12        PC13
## gdp_per_cap                -0.14801344  0.001008133  0.10938078
## inflation                  -0.04385091 -0.004475410  0.05500664
## gni_per_cap                -0.14868511  0.009584271  0.13184416
## urban_p                    -0.05335836 -0.272317139 -0.07970450
## productivity                0.07605659  0.105223475 -0.24213800
## credit_inf                 -0.05331763  0.112913890  0.01080246
## renew_en                   -0.02379598  0.557890932  0.32596461
## renew_electr                0.12373901 -0.200603595 -0.20786441
## electricity_2016            0.14508275  0.072899056  0.17486561
## clean_fuels_2016            0.15368369 -0.057071427  0.11452099
## inf_mort_rate               0.30529794 -0.132938747  0.07541762
## sanitation_services         0.10723198  0.167175678 -0.12344804
## deaths_household_pollution -0.02666143  0.217827160 -0.23058604
## water                       0.26135638  0.117241028  0.48049265
## fert_rate                   0.03304768  0.135884645  0.35453714
## life_exp                   -0.35559877  0.003318500 -0.10672764
## urban_pop_growth            0.27424400  0.137588291 -0.08835846
## emp_agr                    -0.08244732 -0.368030295  0.05243925
## eco_agr                     0.38380045  0.046628532 -0.27886733
## vulner_empl                -0.23914122 -0.405963699  0.39493736
## acc_own                     0.54247332 -0.308554134  0.15171033
##                                    PC14         PC15         PC16
## gdp_per_cap                 0.010711120  0.028195710  0.225007438
## inflation                   0.000495221 -0.033811711 -0.004010458
## gni_per_cap                 0.048615312  0.006540802  0.321001578
## urban_p                    -0.025124159 -0.067700225 -0.076870334
## productivity               -0.176909370  0.083007258 -0.700500594
## credit_inf                  0.086189948 -0.007615905  0.047727022
## renew_en                   -0.354097808  0.047413695 -0.086186809
## renew_electr                0.196196388 -0.042391119 -0.027516007
## electricity_2016           -0.324903114 -0.567810451  0.028543704
## clean_fuels_2016           -0.451291339 -0.068645203  0.056849669
## inf_mort_rate              -0.093309449  0.133294401  0.283603157
## sanitation_services         0.210023728  0.200773886 -0.060755524
## deaths_household_pollution -0.028975365 -0.238341946  0.039196148
## water                       0.291182873  0.406418265 -0.121903272
## fert_rate                   0.484057979 -0.466722652 -0.299024950
## life_exp                    0.173032312 -0.248754800  0.060906438
## urban_pop_growth           -0.067082657  0.103747810  0.188137640
## emp_agr                    -0.242219239  0.131110940 -0.311181149
## eco_agr                     0.119556888 -0.119734612  0.042662069
## vulner_empl                -0.055317974  0.007330463 -0.085152380
## acc_own                     0.042899877 -0.230669674 -0.015777446
##                                    PC17         PC18         PC19
## gdp_per_cap                 0.002153638  0.090949157  0.097326245
## inflation                   0.034938828  0.049119963 -0.007421070
## gni_per_cap                -0.044412343  0.139996813  0.069551130
## urban_p                    -0.150470506  0.105819691 -0.023794930
## productivity                0.184404861 -0.239139516 -0.100347108
## credit_inf                  0.055718134 -0.008838201  0.003339975
## renew_en                   -0.143676048 -0.076769373 -0.171241215
## renew_electr                0.111091441 -0.033903849  0.077475221
## electricity_2016            0.098332618  0.118834906 -0.269289191
## clean_fuels_2016           -0.145724401 -0.090989711  0.520256080
## inf_mort_rate              -0.089926692 -0.133066435 -0.528564451
## sanitation_services        -0.796782837 -0.147643684 -0.093193205
## deaths_household_pollution -0.205124823  0.072134384  0.371511061
## water                       0.246544325  0.175067389  0.092675216
## fert_rate                  -0.146530578  0.100934166  0.085406315
## life_exp                   -0.005764918 -0.111388624 -0.344949381
## urban_pop_growth            0.018623575 -0.003567013  0.062840838
## emp_agr                    -0.254572925  0.625711795 -0.103143970
## eco_agr                     0.118648953 -0.076920538  0.067160273
## vulner_empl                -0.078661483 -0.610078893  0.125719345
## acc_own                    -0.142927505 -0.064926494 -0.024837813
##                                    PC20          PC21
## gdp_per_cap                -0.138153531 -0.7150331379
## inflation                  -0.001232173 -0.0006412151
## gni_per_cap                -0.068321328  0.6907793980
## urban_p                    -0.085943661 -0.0078154078
## productivity                0.066614199  0.0741340121
## credit_inf                  0.024821085 -0.0197353848
## renew_en                    0.007798176 -0.0069109998
## renew_electr                0.016436726  0.0122815014
## electricity_2016           -0.321254707  0.0125601908
## clean_fuels_2016            0.447103471 -0.0071226079
## inf_mort_rate               0.350833901 -0.0340143911
## sanitation_services        -0.187958094 -0.0070090088
## deaths_household_pollution  0.228454492 -0.0114428416
## water                       0.185898660 -0.0045036233
## fert_rate                   0.117816471 -0.0151802643
## life_exp                    0.604475818 -0.0461375741
## urban_pop_growth            0.022260596 -0.0032954164
## emp_agr                     0.109466182 -0.0053926851
## eco_agr                    -0.069884202  0.0032495515
## vulner_empl                -0.132770875  0.0337197537
## acc_own                     0.030706564 -0.0171334248
summary(pca)
## Importance of components:
##                           PC1     PC2    PC3     PC4     PC5     PC6
## Standard deviation     3.5976 1.37042 1.1466 0.96254 0.84946 0.73967
## Proportion of Variance 0.6163 0.08943 0.0626 0.04412 0.03436 0.02605
## Cumulative Proportion  0.6163 0.70576 0.7684 0.81248 0.84684 0.87290
##                           PC7     PC8     PC9   PC10    PC11    PC12
## Standard deviation     0.6735 0.64933 0.55332 0.5422 0.46175 0.44024
## Proportion of Variance 0.0216 0.02008 0.01458 0.0140 0.01015 0.00923
## Cumulative Proportion  0.8945 0.91458 0.92916 0.9432 0.95331 0.96254
##                           PC13    PC14   PC15    PC16   PC17    PC18
## Standard deviation     0.37832 0.37058 0.3549 0.32479 0.3040 0.27949
## Proportion of Variance 0.00682 0.00654 0.0060 0.00502 0.0044 0.00372
## Cumulative Proportion  0.96935 0.97589 0.9819 0.98691 0.9913 0.99503
##                           PC19    PC20    PC21
## Standard deviation     0.23136 0.20584 0.09165
## Proportion of Variance 0.00255 0.00202 0.00040
## Cumulative Proportion  0.99758 0.99960 1.00000

En primer lugar, utilizamos PCA para ver cómo se comportan en alta dimensión nuestras variables. Vemos que el primer componente es muy fuerte, y consigue captar la mayor parte de la varianza total de la muestra, en concreto 3.59 es la desviación estándar de nuestro primer componente. La del segundo es 1.39, y la del tercero 1.13.Los que vienen después captan cada vez menos varianza. Aunque el primero tiene una gran diferencia con respecto a los otros dos (proporción acumulada de varianza en este primer componente es del 61.35%), podría ser buena idea coger 3, ascendiendo hasta el 76.5% de proporción de varianza acumulada. Esto sugiere que hay muchas relaciones ocultas entre nuestras variables y que éstas se pueden explicar en gran medida mediante el uso de 1 componente, aunque con 3 podemos explicar casi un 80% de la varianza total de la muestra.

dfplot <- data.frame(country = countries[-150], region = regions[-150], pc1 = pca$x[ , 1], pc2 = pca$x[ , 2])


qplot(dfplot$pc1,dfplot$pc2,col=dfplot$region,label=dfplot$country,size=I(.1)) + 
  labs(title="First two principal components", x="PC1", y="PC2",color="Region") +
  theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)

En el gráfico superior podemos ver el primer contra el segundo componente en los ejes x e y, con los países y sus regiones que definen el color de la etiqueta. Vemos que claramente países como Luxemburgo, Noruega, Suiza, Islandia o Singapur tienen mucha fuerza en ambos componentes, y son los países que intuitivamente tienen un mayor nivel de vida y por lo tanto esperamos que estén en ese grupo de bien posicionados con respecto a la Second Machine Age y por tanto se estén comiendo una mayor parte de la tarta de los beneficios que el progreso tecnológico provee (bounty).

Como veremos a continuación, los países que aparecen en la parte más alta (alto componente 2) son los más potentes a nivel económico, mientras que los que aparecen más a la derecha alto componente 1) son los mejores en general, los que mejor posicionados están respecto a la escala de spread.

Vamos a dibujar los dos primeros componentes.

barplot(pca$rotation[,1], las=2, col="darkgreen", main = "pc1")

barplot(pca$rotation[,2], las=2, col="darkgreen", main = "pc2")

barplot(pca$rotation[ , 3], las = 2, col = "darkgreen", main = "pc3")

El primer componente vemos que utiliza casi todas las variables, consigue explicar una parte de todas. Tiene peso positivo para el gdp per capita, el gni per capita, urban population, productivity, credit information, electricity 2016, clean fuels 2016, sanitation services, water, life expectancy y account ownership. Por otro lado, tiene pesos negativos para inflation, renewable energy, renewable electricity, infant mortality rate, deaths for household pollution, fertility rate, urban population growth, emp agr (empleo en agricultura – porcentaje del labour force dedicada a la agricultura), economy agriculture (porcentaje de la economía doméstica proveniente de la agricultura) y vulnerable employment. Por lo tanto podríamos decir que este primer componente es algo parecido a lo que venimos buscando, pues explica bastante bien las diferencias entre los países, desde económicas hasta sanitarias pasando por el tipo de energía que utilizan.

El segundo componente tiene un peso alto positivo para gdp per capita, gni per capita, urban population, productivity, energías renovables, electricidad renovable, negativo para electricity 2016, negativo para clean fuels 2016, negativo para el agua (muy bajo; no le da mucho peso ni al acceso al agua ni a los servicios sanitarios), positivo para urban pop growth y positivo para account ownership. Vemos por tanto que este segundo componente es principalmente un componente económico puramente. En él España por ejemplo tiene menos fuerza que otros países como Irlanda (aunque siendo un paraíso fiscal tiene trampa esta conclusión).

Por último, el gráfico 3 nos muestra que el tercer componente contiene principalmente las variables credit information, renewable energy y renewable electricity, es decir, cómo de avanzado está el país en cuanto a sostenibilidad y el nivel de conocimiento general que hay del crédito, lo cual es fundamental para que un país aproveche la ola del progreso tecnológico, pues sin una participación activa del sistema económico, un menor porcentaje de la población tendrá acceso a los beneficios de ese desarrollo. El uso de energías limpias y la cultura general de la población son fundamentales para que el crecimiento económico y ese race against the machine del país sea sostenible.

Después de este PCA, vamos a utilizar FA para finalizar con el estudio y definir formalmente la variable spread, la cual posiblemente encontremos contenida entre varios factores distintos. Vamos a incluir 3 factores, pues como hemos visto arriba el uso de 3 componentes principales nos permitiría explicar el 76% de la varianza de la muestra y, lo que es mejor, los 3 son interpretables fácilmente y es muy sencillo aplicarlos a la cuestión que tratamos de resolver. Viendo el screeplot bien podríamos haber seguido el criterio del codo y haber cortado en 1, pero va a ser insuficiente para explicar las diferencias. Esperamos que el primer factor se encargue de explicar el spread de forma general, atendiendo a todas las variables que esperamos sean explicadas por esta variable latente; los otros dos factores posiblemente, basándonos en los resultados de PCA (si utilizamos FA sin rotación) sean el componente del spread que tiene más que ver con lo puramente económico y el componente del mismo que tiene más que ver con cuestiones de sostenibilidad y participación ciudadana del sistema.

3.2. Factor Analysis

No aplicaremos ninguna rotación, pues tras haber probado con varios tipos nos encontramos con que rotation = “none” es el método que nos da unos factores más interpretables, además de acumular la mayor parte de la varianza captada por los factores en el primer factor, lo cual es positivo para nuestro objetivo principal que es encontrar la variable latente spread.

factors <- factanal(datos_brutos, factors = 3, rotation = "none", scores = "regression", lower = 0.1)

factors
## 
## Call:
## factanal(x = datos_brutos, factors = 3, scores = "regression",     rotation = "none", lower = 0.1)
## 
## Uniquenesses:
##                gdp_per_cap                  inflation 
##                      0.100                      0.846 
##                gni_per_cap                    urban_p 
##                      0.100                      0.443 
##               productivity                 credit_inf 
##                      0.123                      0.664 
##                   renew_en               renew_electr 
##                      0.316                      0.811 
##           electricity_2016           clean_fuels_2016 
##                      0.164                      0.130 
##              inf_mort_rate        sanitation_services 
##                      0.111                      0.114 
## deaths_household_pollution                      water 
##                      0.252                      0.188 
##                  fert_rate                   life_exp 
##                      0.177                      0.100 
##           urban_pop_growth                    emp_agr 
##                      0.526                      0.114 
##                    eco_agr                vulner_empl 
##                      0.311                      0.129 
##                    acc_own 
##                      0.320 
## 
## Loadings:
##                            Factor1 Factor2 Factor3
## gdp_per_cap                 0.736   0.635         
## inflation                  -0.367          -0.138 
## gni_per_cap                 0.745   0.622         
## urban_p                     0.719   0.131  -0.152 
## productivity                0.816   0.444  -0.116 
## credit_inf                  0.520  -0.169   0.192 
## renew_en                   -0.719   0.318   0.257 
## renew_electr               -0.224   0.173   0.330 
## electricity_2016            0.839  -0.361         
## clean_fuels_2016            0.899  -0.212  -0.130 
## inf_mort_rate              -0.889   0.198  -0.242 
## sanitation_services         0.907  -0.250         
## deaths_household_pollution -0.853   0.117         
## water                       0.857  -0.267         
## fert_rate                  -0.850   0.239  -0.208 
## life_exp                    0.917           0.254 
## urban_pop_growth           -0.634   0.269         
## emp_agr                    -0.895   0.101   0.272 
## eco_agr                    -0.801           0.197 
## vulner_empl                -0.903           0.225 
## acc_own                     0.788   0.219   0.101 
## 
##                Factor1 Factor2 Factor3
## SS loadings     12.683   1.732   0.643
## Proportion Var   0.604   0.082   0.031
## Cumulative Var   0.604   0.686   0.717
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 768.77 on 150 degrees of freedom.
## The p-value is 8.13e-84
cbind(factors$loadings, factors$uniquenesses)
##                               Factor1     Factor2       Factor3          
## gdp_per_cap                 0.7359662  0.63452938  4.155716e-02 0.1000000
## inflation                  -0.3665357 -0.02282047 -1.383439e-01 0.8459872
## gni_per_cap                 0.7445106  0.62211085  4.887288e-02 0.1000000
## urban_p                     0.7187899  0.13091917 -1.523659e-01 0.4429849
## productivity                0.8161107  0.44424276 -1.156209e-01 0.1232401
## credit_inf                  0.5196299 -0.16925103  1.921919e-01 0.6644133
## renew_en                   -0.7188026  0.31770524  2.569784e-01 0.3163405
## renew_electr               -0.2236911  0.17329489  3.299529e-01 0.8111107
## electricity_2016            0.8392060 -0.36107547  4.217055e-02 0.1635819
## clean_fuels_2016            0.8989310 -0.21237754 -1.295020e-01 0.1300528
## inf_mort_rate              -0.8894618  0.19823458 -2.421372e-01 0.1109326
## sanitation_services         0.9071072 -0.24988358  1.530662e-02 0.1144710
## deaths_household_pollution -0.8526162  0.11689647  8.460251e-02 0.2522163
## water                       0.8569504 -0.26721429  7.670451e-02 0.1883512
## fert_rate                  -0.8501949  0.23862068 -2.082418e-01 0.1768658
## life_exp                    0.9171641 -0.06917325  2.536867e-01 0.1000000
## urban_pop_growth           -0.6335769  0.26899404 -2.472499e-05 0.5262214
## emp_agr                    -0.8951666  0.10121577  2.724716e-01 0.1141917
## eco_agr                    -0.8005125  0.09459981  1.968513e-01 0.3114842
## vulner_empl                -0.9034610  0.06565532  2.245987e-01 0.1290051
## acc_own                     0.7883023  0.21906972  1.007613e-01 0.3204361

Vemos que con estos 3 factores podemos explicar el 71.6% de la varianza total de la muestra, siendo la proporción de varianza que explica cada factor la siguiente:

Factor 1: 60.3%

Factor 2: 8.2%

Factor 3: 3.1%

Viendo las uniquesses de cada variable, la mayoría son bastante bajas, lo cual nos indica que estos 3 factores explican la mayoría de esas variables. Sin embargo, hay algunas más altas como las de inflation (0.84), urban population (0.44), credit information (0.657), renewable electricity (0.81), urban population growth (0.52). Para mejorar el modelo, podríamos eliminar estas variables, pues están introduciendo ruido a la hora de encontrar la variable latente spread.

barplot(factors$loadings[,1], las = 2, col = "darkgreen", main = "factor1")

Vemos que el primer factor sería, como esperábamos, uno que tiene en cuenta todas las variables prácticamente. Podemos identificar este primer factor como el spread, la variable principal que estábamos buscando.

barplot(factors$loadings[,2], las = 2, col = "darkgreen", main = "factor1")

barplot(factors$loadings[,3], las = 2, col = "darkgreen", main = "factor1")

El segundo factor sería el componente económico del spread, es decir cuál es la parte económica del “bounty” (beneficios debidos al progreso tecnológico), mientras que el tercer factor sería el componente calidad de vida más allá de lo económico, pues como vemos tiene un peso fuerte para las variables credit information, renewable energy, renewable electricity, negativo para infant mortality rate, fertility rate; y por último positivo para emp agr, life expectancy, eco agr, vulnerable employment y account ownership. Lo que menos nos cuadraría sería el vulnerable employment teniendo una relación positiva con el factor calidad de vida más allá de lo económico (o sostenibilidad como dijimos antes), pues esperaríamos que tuvieran una relación negativa.

dfplot <- data.frame(country = countries[-150], region = regions[-150], factor1 = factors$scores[ , 1], factor2 = factors$scores[ , 2])


qplot(dfplot$factor1,dfplot$factor2,col=dfplot$region,label=dfplot$country,size=I(.1)) + 
  labs(title="First two factors", x="SPREAD IN GENERAL (FACTOR 1)", y="ECONOMIC PART OF SPREAD (FACTOR 2)",color="Region") +
  theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)

Es especialmente llamativo el caso de Luxemburgo, que es claramente la que mejor posicionada está en el spread; tanto en el general como en la parte más económica del mismo. Los países que están más a la derecha y más arriba son los que se están repartiendo el pastel del progreso tecnológico principalmente. Vemos que coincide con los resultados que podríamos esoerar de forma intuitiva, con países muy desarrollados marcando el ritmo de esta carrera por el progreso: Singapur, Noruega, Islandia… Concluimos por tanto que hemos encontrado lo que veníamos buscando, un factor que pudiéramos identificar como la variable spread. La hemos encontrado realmente repartida en 3 factores: spread general, spread tecnológico y spread sostenible.

Probaremos, por último, a realizar el mismo análisis eliminando las variables que tengan uniquesses demasiado altas, para tratar de afinar un poco más la puntería. De momento nos quedaríamos con el modelo presentado hasta ahora.

datos_reduced <- data.frame(datos_brutos)

datos_reduced[ , c("inflation", "urban_p", "credit_inf", "renew_electr", "urban_pop_growth")] <- NULL

factors_reduced <- factanal(datos_reduced, factors = 3, rotation = "none", scores = "regression", lower = 0.1)

factors_reduced
## 
## Call:
## factanal(x = datos_reduced, factors = 3, scores = "regression",     rotation = "none", lower = 0.1)
## 
## Uniquenesses:
##                gdp_per_cap                gni_per_cap 
##                      0.100                      0.100 
##               productivity                   renew_en 
##                      0.126                      0.337 
##           electricity_2016           clean_fuels_2016 
##                      0.162                      0.122 
##              inf_mort_rate        sanitation_services 
##                      0.107                      0.112 
## deaths_household_pollution                      water 
##                      0.244                      0.191 
##                  fert_rate                   life_exp 
##                      0.188                      0.100 
##                    emp_agr                    eco_agr 
##                      0.129                      0.307 
##                vulner_empl                    acc_own 
##                      0.127                      0.322 
## 
## Loadings:
##                            Factor1 Factor2 Factor3
## gdp_per_cap                 0.738   0.632         
## gni_per_cap                 0.747   0.620         
## productivity                0.817   0.441  -0.114 
## renew_en                   -0.716   0.316   0.224 
## electricity_2016            0.839  -0.365         
## clean_fuels_2016            0.899  -0.218  -0.149 
## inf_mort_rate              -0.891   0.202  -0.244 
## sanitation_services         0.908  -0.253         
## deaths_household_pollution -0.855   0.125         
## water                       0.854  -0.269         
## fert_rate                  -0.848   0.235  -0.193 
## life_exp                    0.918           0.263 
## emp_agr                    -0.891   0.102   0.257 
## eco_agr                    -0.799           0.213 
## vulner_empl                -0.902           0.234 
## acc_own                     0.788   0.223         
## 
##                Factor1 Factor2 Factor3
## SS loadings     11.302   1.586   0.445
## Proportion Var   0.706   0.099   0.028
## Cumulative Var   0.706   0.806   0.833
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 536.32 on 75 degrees of freedom.
## The p-value is 7.69e-71
cbind(factors_reduced$loadings, factors_reduced$uniquesses)
##                               Factor1     Factor2       Factor3
## gdp_per_cap                 0.7383910  0.63207775  0.0363344297
## gni_per_cap                 0.7466168  0.61993000  0.0434230561
## productivity                0.8168544  0.44062459 -0.1138407175
## renew_en                   -0.7164421  0.31627830  0.2237601169
## electricity_2016            0.8387445 -0.36521045  0.0306045157
## clean_fuels_2016            0.8989023 -0.21755158 -0.1488108709
## inf_mort_rate              -0.8906982  0.20172388 -0.2437547718
## sanitation_services         0.9077532 -0.25303007  0.0003041284
## deaths_household_pollution -0.8545754  0.12484377  0.0988855429
## water                       0.8544821 -0.26897616  0.0796860173
## fert_rate                  -0.8481071  0.23528000 -0.1933136076
## life_exp                    0.9182668 -0.07396433  0.2625352882
## emp_agr                    -0.8911070  0.10239945  0.2569312519
## eco_agr                    -0.7988382  0.09736057  0.2126371243
## vulner_empl                -0.9022263  0.06629560  0.2344317865
## acc_own                     0.7878089  0.22293001  0.0884406168

Vemos que el modelo funciona bastante mejor, pues los 3 factores conjuntamente explican un 83.1% de la varianza de la muestra, mientras que el primer factor explica el 70.5%. Como dijimos antes, el p-value se usa en la práctica para comparar modelos. Si comparamos el p-value en cada modelo, vemos que en este segundo modelo después de eliminar esas variables que estaban introduciendo ruido en la muestra el p-value baja significativamente. Podríamos quedarnos, pues, con este modelo, pues descartamos que los factores que componen el spread (spread general + spread económico + spread sostenible) sean variables latentes que expliquen la inflación y el resto de variables que hemos eliminado.

Dibujamos los gráficos correspondientes a este último modelo:

dfplot <- data.frame(country = countries[-150], region = regions[-150], factor1 = factors_reduced$scores[ , 1], factor2 = factors_reduced$scores[ , 2])


qplot(dfplot$factor1,dfplot$factor2,col=dfplot$region,label=dfplot$country,size=I(.1)) + 
  labs(title="First two factors", x="SPREAD IN GENERAL (FACTOR 1)", y="ECONOMIC PART OF SPREAD (FACTOR 2)",color="Region") +
  theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)

Vemos que Chad, Somalia o Sierra Leona son algunos de los países que tienen en el spread general una puntuación muy negativa, sin embargo en el spread económico en concreto no salen tan mal parados, y tienen una puntuación similar a países como Australia.

barplot(factors_reduced$loadings[,1], las = 2, col = "darkgreen", main = "factor1")

barplot(factors_reduced$loadings[,2], las = 2, col = "darkgreen", main = "factor2")

barplot(factors_reduced$loadings[,3], las = 2, col = "darkgreen", main = "factor3")

Vemos que los factores que nos salen, por la matriz L que estamos dibujando vector a vector arriba (del factor 1 al factor 3), son muy similares a los obtenidos cuando incluíamos todas las variables. La mejora viene del lado de la reducción del ruido. Al trabajar con 16 variables en lugar de con 21 podemos extraer más información de la muestra y así definir más claramente la variable spread.

4. Visualización de Resultados

A continuación dibujaremos los resultados obtenidos sobre un mapa del mundo, mostrando el score de cada país en cada uno de los 3 factores, que recordemos son:

Primero dibujamos los gráficos para el primer modelo:

if(!require(rworldmap)) {
  
  install.packages("rworldmap")
  
}

matched <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,1]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
matched2 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,2]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
matched3 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,3]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
mapCountryData(matched, nameColumnToPlot="value", mapTitle="General Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")

mapCountryData(matched2, nameColumnToPlot="value", mapTitle="Economic Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")

mapCountryData(matched3, nameColumnToPlot="value", mapTitle="Sustainable Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")

# zoom in Africa
mapCountryData(matched, nameColumnToPlot="value", mapRegion = "Africa", mapTitle="General Spread by Factor Analysis: Zoom in Africa", 
               catMethod = "pretty", colourPalette = "topo")

#zoom in Europe
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Europe", mapTitle = "General Spread by Factor Analysis: Zoom in Europe", catMethod="pretty", colourPalette = "topo")

#zoom in America
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "North America", mapTitle = "General Spread by Factor Analysis: Zoom in North America", catMethod="pretty", colourPalette = "topo")

mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Latin America", mapTitle = "General Spread by Factor Analysis: Zoom in South America", catMethod="pretty", colourPalette = "topo")

#zoom in Asia
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Asia", mapTitle = "General Spread by Factor Analysis: Zoom in Asia", catMethod="pretty", colourPalette = "topo")

Vemos que el factor verdaderamente revelador es el del general spread; es el que mejor nos consigue explicar las diferencias entre países y es por lo tanto el resultado al que más atención deberemos prestar en nuestro estudio. Podemos ver que el mundo occidental en general, con América del Norte, Australia, Europa (Western), y algunas regiones aisladas (Japón, United Arab Emirates), son los países que tienen un mejor posicionamiento con respecto a los cambios que está sufriendo el mundo debido al progreso tecnológico, y por lo tanto son los países que más valor están absorbiendo del que se está creando (bounty). Vemos que en definitiva estas dos variables (spread y bounty) están muy relacionadas; el spread es el que mide cómo se distribuye el bounty en el mundo, es decir es una escala que nos permite determinar qué países están aventajando a los demás debido al progreso, qué países están bien posicionados en la balanza cada vez menos equilibrada de la calidad de vida relativa en el mundo. Vemos que existe una desigualdad estructural demográfica, lo cual nos viene a decir que los países colonialistas del siglo pasado son también, en su mayoría, los que están liderando el cambio tecnológico, por lo tanto distanciándose más del resto de países. Especialmente en África observamos que es un continente que se está quedando atrasado respecto al progreso, ya que muchos países de este contienente tienen un score verdaderamente bajo en el spread general. Este mapa nos sirve también para relacionar esta desigualdad derivada del progreso tecnológico con la situación política internacional. Vemos que algunos países como Israel, pese a estar en una zona en la que la mayoría de los países tienen un score mediano aproximadamente en el factor 1, tienen un score equiparable al de Europa, EEUU o Australia, debido a sus conexiones políticas internacionales.

Los mejores países para vivir en líneas generales serían Irlanda (que se cuela entre los países con un score más alto en spread general), Noruega, Suiza, Luxemburgo, Singapur… Estos son los países que mejor puntúan en el spread general, y por lo tanto los que más se están distanciando del resto.

En cuanto al spread económico puro, podemos ver que destacan muy pocos países en este aspecto. Estos son Suiza, Noruega, Islandia e Irlanda. De forma negativa vemos a Egipto y a otros países del norte de África, y algunos países de Europa del Este.

Respecto al Spread Sostenible, vemos que el país que mejor puntúa sobre el resto de forma notable es Nepal, lo cual encaja debido a su religión; precisamente el budismo se centra en buscar el equilibrio; por lo que tiene sentido que también crezcan de forma sostenible (con pocas muertes infantiles, acceso al agua y a la sanidad, alta esperanza de vida, alto uso de energías renovables…). Además, tener un porcentaje importante de economía agrícola podría verse como una forma de crecer más sostenible, pues es una economía enfocada en el alimento, además de acarrear generalmente una mayor calidad de la comida doméstica y ser ésta más saludable.

Dibujaremos por último los mismos gráficos para el segundo modelo de factores, aunque esperamos que los resultados sean prácticamente los mismos como se explica en la anterior sección, pues cambia muy poco la interpretación de los mismos.

matched <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,1]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
matched2 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,2]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
matched3 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,3]), joinCode="NAME", nameJoinColumn="country")
## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data
mapCountryData(matched, nameColumnToPlot="value", mapTitle="General Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification

mapCountryData(matched2, nameColumnToPlot="value", mapTitle="Economic Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")

mapCountryData(matched3, nameColumnToPlot="value", mapTitle="Sustainable Spread by Factor Analysis", 
               catMethod = "pretty", colourPalette = "topo")

# zoom in Africa
mapCountryData(matched, nameColumnToPlot="value", mapRegion = "Africa", mapTitle="General Spread by Factor Analysis: Zoom in Africa", 
               catMethod = "pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification

#zoom in Europe
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Europe", mapTitle = "General Spread by Factor Analysis: Zoom in Europe", catMethod="pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification

#zoom in America
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "North America", mapTitle = "General Spread by Factor Analysis: Zoom in North America", catMethod="pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification

mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Latin America", mapTitle = "General Spread by Factor Analysis: Zoom in South America", catMethod="pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification

#zoom in Asia
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Asia", mapTitle = "General Spread by Factor Analysis: Zoom in Asia", catMethod="pretty", colourPalette = "topo")
## You asked for 7 categories, 10 were used due to pretty() classification